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ABSTRACT 

Galaxies with stellar bulges are generically observed to host supermassive black holes 
(SMBHs). The hierarchical merging of galaxies should therefore lead to the formation 
of SMBH binaries. Merging of old massive galaxies with little gas promotes the forma- 
tion of low-density nuclei where SMBH binaries are expected to survive over long times. 
If the binary lifetime exceeds the typical time between mergers, then triple black hole 
(BH) systems may form. We study the statistics of close triple-SMBH encounters in 
galactic nuclei by computing a series of 3-body orbits with physically-motivated initial 
conditions appropriate for giant elliptical galaxies. Our simulations include a smooth 
background potential consisting of a stellar bulge plus a dark matter halo, drag forces 
due to gravitational radiation and dynamical friction on the stars and dark matter, 
and a simple model of the time evolution of the inner density profile under heating 
and mass ejection by the SMBHs. We find that the binary pair coalesces as a result of 
repeated close encounters in ~85% of our runs, and in ~15% of cases a new eccentric 
binary forms from the third SMBH and binary remnant and coalesces during the run 
time. In about 40% of the runs the lightest BH is left wandering through the galactic 
halo or escapes the galaxy altogether, but escape of all three SMBHs is exceedingly 
rare. The triple systems typically scour out cores with mass deficits ~l-2 x their total 
mass, which can help to account for the large cores observed in some massive elliptical 
galaxies, such as M87. The high coalescence rate, prevalence of very high-eccentricity 
orbits, and gravitational radiation "spikes" during close encounters in our runs, may 
provide interesting signals for the future Laser Interferometer Space Antenna (LISA) . 

Key words: black hole physics — cosmology: theory — galaxies: elliptical and lentic- 
ular, cD — galaxies: interactions — galaxies: nuclei — methods: numerical 



1 INTRODUCTION 

In the favored cold dark matter cosmology, present-day 
galaxies were assembled hierarchically from smaller build- 
ing blocks at earlier cosmic times. Since all nearby galax- 
ies with stellar spheroids are observ ed to host nuclear 
SMBHs (|Kormendv fc Gebhardtl 1200 ll h hierarchical merg- 
ing leads inevitably t o the formation of SMBH binaries 
ijBegelman et alj Il980t ) . If the binary lifetime exceeds the 
typical time between mergers, then some galactic nuclei 
should contain systems of three or more SMBHs. These sys- 
tems are particularly interesting as they often lead to the 
ejection of one of the B Hs at a speed compara ble to the 
galactic escape velocity (|Hoffman fc Loeq l2006h . In mas- 
sive elliptical galaxies the typical speeds are ~ 10 3 km s — 1 , 
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far gre ater than attainable through gravitational radiation 
recoil (jCentrellal 120061 : iFavata et aiTl2004l | Blanchet et al.l 




Spatially resolved pairs of nuclei have been observed in a 
few active galaxies. The most famous example is NGC 6240, 
an Ultraluminous Infrared Galaxy (ULIRG) in which two 
distinct active galactic nuclei (AGN) are clea rly seen in hard 
X-ray s at a project e d separation of ~1 kpc jKomossa et al.l 
l2003n . iMaoz et al.1 l| 19951 . l2005h observed a variable UV 
source, possibly a second active nucleus, at a projected sep- 
aration of ~ 60 pc from the primary nucleus in the spiral 
galaxy NGC 47 3 6, wh ich shows signs of a recent merger. 
iRodriguez et all (|2006l ) have detected what is thought to be 
an SMBH binary at a projected separation of just 7.3 pc in 
the radio galaxy 0402+379, through multi-frequency radio 
observations using the Very Long Baseline Array (VLB A). 
We begin by discussing the theory of how such systems 
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evolve, and the conditions under which they might acquire 
a third BH. 



1.1 Black hole binaries 

When two galaxies merge, their dense nuclei sink to the 
center of the merger product by dynamical friction. As the 
nuclei spiral in, tidal forces gradually strip the two SMBHs 
of their surrounding stars and dark matter. In mergers be- 
tween galaxies of comparable mass, the BHs are able to come 
together and form a bound SMBH binary on a timescale of 
order 10 9 yrs. The binary continues to harden by dynamical 
friction until it reaches a separation of order 
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known as the "hardening radius" (e.g. lOurn lan 1996). Here 
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(i = m 1 m2/{m\ + m-i) is the reduced mass of the two BHs 
with masses mi and mi, a is the velocity dispersion of the 
stars beyond the binary's sphere of gravitational influence, 
q is the binary mass ratio mi/mi 1, and m oin — m 1 + mi 
is the total mass of the binary. For smaller separations the 
binary looks like a point mass to the distant stars contribut- 
ing to dynamical friction, but close stellar encounters pref- 
erentially harden the binary and so dominate further energy 
loss. Only stars on nearly radial orbits, with periapsis dis- 
tances of order the binary separation, can extract energy 
from ("harden") the binary in this stage. These stars un- 
dergo strong 3-body interactions with the binary and es- 
cape its vicinity with speeds comparable to the black holes' 
orbital speed. In the low-density nuclei of large elliptical 
galaxies, the total mass in stars on such "loss cone" orbits 
is small compared to the mass of the binary. Furthermore 
the two-body stellar relaxation time is long compared to a 
Hubble time, so once the stars initially on l oss cone orbits 
are c l eared out, the loss cone rem ai ns empty (|Frank fc Reesl 
1 19761 ; iLightman fc Shapird 1 19771 : I Cohn fc Kulsrudl 1 19781 ). 
Since the binary must eject of order its own mass per e- 
folding in its semi-major axis, the system stops hardening 
around ahard unless some other mechanism causes sufficient 
mass flux through the binary. 

If the binary reaches a separation around 
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10 10 yrs, 

where e is the orbital eccentricity of the binary and /(e) = 
(1 - e 2 ) 7/2 /(l + 73e 2 /24 + 37e 4 /96), then it can coa- 
lesce on a timescale T aw through gravitational radiation 
l|Begelman et al.lll980l ). To get from a hard to a gw it must 
bridge a gap 
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by some mechanism other than stellar-dynamical friction or 
gravitational radiation. The question of whether and how 
it crosses this gap has become known a s the "final parsec 
problem" (|Merritt fc Milosavlievil 120051 ). 

In many galaxies there probably are alternative mecha- 
nisms for crossing the gap. When gas-rich galaxies merge, 
tidal torques channel large amounts of gas into the cen- 



tral ~ 100 pc jBvrd et al.lll987l ; lHernquistlll989l ). The gas 
may lose energy through radiation and angular momen- 
tum through viscous torques, and is therefore not subject 
to a loss cone proble m. Using Smoothed Part icle Hydro- 
dynamics simulations lEscala et al.l (|2004l . l2005h compute a 
merger time of order 10 7 yrs in an environment typical of 
the central regions of ULIRGs, which are tho ught to be gas- 
rich galaxies caught in the act of merging (jSanders et al.l 
1988). The nuclei of galaxies are also observed to con- 
tain numerous massive perturbers (MPs) such as star clus- 
ters, molecular clouds, and possibly intermediate- mass black 
holes (IMBHs) . These objects scatter stars into the loss cone 
much more efficiently than other stellar mass objects, since 
the relaxation rate scales as the perturber mass for a fixed 
mass density of perturbers. IPerets et al.l (|2006 | ) extended 
the Fokker-Planck loss cone formalism dFrank fc Reeslll976l ; 
ILightman fc Shapird [l977l ; ICohn fc Kulsrudl 1 19781 ) to acco- 
modate a spectrum of perturber masses and account for re- 
laxation by rare close encounters with MPs. They show that 
the population of known MPs in the nucleus of the Milky 
Way is sufficient to bring a 4 x 10 6 Mq BH binary to a gw in 
~ 6x 10 s yrs, and it is reasonable to expect similar perturber 
populations in other star-forming spiral galaxies. 

The final parsec problem is often mentioned as a 
caveat when predicting the SMBH coalescence signal in low- 
frequency gravitational wave detectors such as the upcoming 
Laser Interferometer Space Anntena (LISA). However the 
LISA event rate is exp ected to be dominated by small galax- 
ies a t high redshift (IWvithe fc Loebl l2003al ; ISesana et al.l 
120051 : iRhook fc Wvithell2005h . where the gas content and 
central densities tend to be high and the relaxation times 
short. For this reason the stalling problem is probably not a 
significant concern for the LISA SMBH coalescence signal. 
On th e other hand BH ejections by grav itational radiation 
recoil jMerritt et alj|2004 lHaimanll2004l ) may play an im- 
portant role in the high-redshift coalescence rate. The long- 
term survival of SMBH binaries is likewise unlikely in the 
gas-rich cores of quasars and ULIRGs. 

However none of the gap-crossing mechanisms discussed 
so far are likely to reduce the coalescence time below a 
Hubble time in mergers betwe en giant, gas-poor elliptical 
galaxies. iMerritt fc Pood l|2004l ) show that a significant frac- 
tion of stars on "centrophilic" orbits in a triaxial potential 
can greatly increase the mass flux into the loss cone. Some 
non-axisymmetric potentials can also excite bar instabili- 
ties that cause rapi d mass flow throug h the binary and ef- 
ficient coalescence l|Berczik et al.l 12006 ) . However a central 
SMBH can disrupt box orbits and ind uce axisymmetry in 
the inner regions of a triaxial galax y (|Merritt fc Quinlad 
1 199a ; iHollev-Bockelmann et al.l 120021 ). and it is uncertain 
how often these geometry-specific mechanisms bring the co- 
alescence time below a Hubble time. 

One can naively assess the likelihood of coalescence 
by considering the "full" and "empty" loss cone hardening 
times, Tfn,u and T emp ty, in the nuclei of various galaxies as- 
suming a spherical and isotropic distribution function. Tf u u 
is the hardening time assuming every star kicked out of the 
loss cone is instantly replaced, while t< 



empty 



is the time as- 



suming stellar two-body relaxation to be the only replen- 
ishing mechanism. In small, dense galaxies Tf u u ~ 10 5-6 
yrs and 1 empty ~ 10' yrs while in the lowest-density 
cores of giant ellipticals and cD galaxies Tf u u ~ 10 8 yrs 
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and T empt y ~ 10 14 yrs (Yu 2002). While the empty loss cone 
rate is difficult to believe in any galaxy given at least some 
clustering on scales larger than 1 M©, it also seems difficult 
to approach the full loss cone rate if there is no gas around 
and no strong radial bias in the stellar distribution. From 
this point of view the stalling of binaries seems unlikely in 
small galaxies but probable in low-density, gas-poor ellipti- 
cals. If some binaries do survive for around a Hubble time, 
then the hierarchical buildup of galaxies will inevitably place 
three or more SMBHs in some merging systems. 



1.2 Merger-induced binary evolution before 
3-body interactions: Back-of-the-envelope 
calculations 

An inspiralling satellite affects the evolution of a binary 
SMBH even long before it sinks to the center, by perturb- 
ing the large-scale potential and scattering stars into the 
loss cone. We may estimate the extent of this effect as a 
function of satellite mass and distance from the center of 
the host ga laxy using a rough but simple argument due to 
iRoos! l|l98ll) . The change in velocity necessary to deflect a 
star at radius q into the loss cone is AV ~ hi c /q, where 
hi c ~ u^/ri n frbin is the characteristic specific angular mo - 
mentum of stars on loss cone orbits |Frank fc Reeslll976h . 
Tinj = Gmi, irl /a' 2 is the SMBHs' radius of influence, and 
Thin is the binary separation. The dynamical time at this 
radius is tdyn ~ so the acceleration required to scat- 

ter a star into the loss cone is roughly ai c (q) ~ AV/td yn ~ 
v 2 y/rinfTbin/ q 2 ■ Equating this with the tidal acceleration 



caused by the satellite, a t id = 2GM, 
the satellite's radius, yields q 3 = a 2 
or with r bin = a hard = Gp/4a 2 , 



(r)q/r 3 , where r is 
y/ri„fr b i„r 3 /2GM aat (r), 
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L4M sat (r)_ 

The r-dependence of M sa t reflects the tidal stripping of the 
satellite as it spirals inward. Equation @ defines a criti- 
cal radius q, outside of which the satellite can deflect stars 
into (and out of) the loss cone in one dynamical time. The 
mass flux through the binary induced by the satellite is then 
approximately 

dM 3 



dt 
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where p{q) is the density of the host galaxy at radius q and 
0f c ~ ri n frbin/r 2 is the geometrical factor accounting for the 
fraction of stars on loss cone orbits as a functi on of radius r, 
assum ing an isotropic distribution function (|Frank fc Reesl 
1 19761 ). For a fixed satellite mass and distance, we can then 
define a "binary feeding" timescale by 

Tfeed = dM stars /dt = 2tv p(q)q^ae 2 c (q) ' W 

To determine whether the scattering of stars into the 
loss cone by the satellite is sufficient to harden the binary 
enough to prevent a close 3-body encounter before the in- 
truder arrives at the galactic center, we must compare Tf ee d 
with the timescale on which the satellite spirals in by dy- 
namical friction. In the approximation of slow inspiral we 
may write the dynamical friction timescale as Tdf = \r/r\ ~ 
\v/(dv/dt)df I- Substituting (dv/dt)df from Chandrasekhar's 




r (kpc) 

Figure 1. Comparison of the "feeding timescale" , Tf ee d, on which 
an inspiralling satellite scatters mass into the loss cone of an 
SMBH binary, with the dynamical friction timescale, r d f, on 
which the satellite spirals in. Upper (red dashed) line: Tf ee£ j com- 
puted from equation JB); Lower curve: tm computed from equa- 
tion Q. Both timescales are plotted as a function of the satellite's 
distance from the center of the host galaxy. For this plot we chose 
a binary mass of 4.5 X 10 8 M© and merger mass ratio of 3:1 in 
the stars. The galactic model is discussed in the text. 



formula (equation (128 p in §2.4.1; [Chandrasekhar 1943) yields 
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where X = v sa t/V2a. v sa t(r) in equation is com- 
puted from v sa t(r) = \/GM host (r)/r, where M ho st(r) = 
M s tars{r) + Mhalo(r) + m b in is the mass of the host galaxy 
enclosed within radius r. M sa t(r) is the satellite mass con- 
tained within the tidal truncation radius obtained from a 
simple point mass approximation, r t id ~ [Af sat /Mh 03t ] 1//3 r 
(this slightly underestimates r t i d as the satellite approaches 
the center of the host). In Fig. 1 we plot Tdf and Tf ee d 
as a function of r for a satellite with one third the stellar 
mass of the host, which contains a binary with (7711,7712) = 
(1.2,3.7) x 10 s Mq . Both host and satellite are modelled as 
Hernquist profiles (|Hernquistlll990l) . with their masses and 
effective radii set by observed scaling relations. The details 
of the galactic model are described further in §2.1-2. 

Since Tf ee d remains about an order of magnitude above 
Tdf throughout the inspiral, this simple calculation makes it 
plausible that the binary survives the merger process and 
undergoes close triple interactions with the infalling SMBH. 
The tidal approximation (as well as our treatment of dynam- 
ical friction) breaks down as the satellite approaches r in f, 
so the plot is cut off at a separation of ~ 100 pc, when the 
satellite still has ~4 e- foldings to go to reach a^ard- How- 
ever this final stage of the inspi ral is found to proceed very 
rapidly in N-body simulation s dQuinlan fc Hernquistl Il997l ; 
iMilosavlievic fc Merritt|[200ll ; iMerrittl 12006). The merger's 
effect on the binary may be dominated by viol ent relaxation 
or collective effects such as a bar instability (jBerczik et al.l 
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2006), in which case our two-body approach does not cap- 
ture its essence. The evolution of the core distribution func- 
tion under the influence of a major merger is an intriguing 
open problem for simulators. 

After the third BH becomes bound to the binary (but 
still before the onset of close 3-body interactions) another 
hardening mechanism may become important. If the angle 
of inclination i of the outer binary (formed by the intruder 
and the inner binary center-of-mass (COM)) exceeds a crit- 
ical angle 6 cr it ~ 39°, then the quadrupolar perturbation 
from the in truder induc es eccentricity oscillations through a 
maximum |Kozailll962l ) 



(8) 



Since the gravitational radiation rate increases sharply 
toward high eccentricities, these "Kozai oscillations" can 
greatly enhance the radiation, possibly causing the binary 
to coalesce before it can undergo strong 3-body interac- 
tions with the intruder dBlaes et al.l 120021 ) . General rela- 
tivistic precession can d estroy the Kozai resonance (e.g. 
iHolman et ail 1 19971 ) . but iBlaes et~aH l|2002l ) find that this 
does not happen for 
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where a,i n and a ou t are the semi-major axes of the inner and 
outer binaries, q out is the outer binary mass ratio, and e; n 
and e ut are the inner and outer eccentricities. This leaves a 
window of about a factor of 10 in a out /a,i n in which the Kozai 
mechanism can operate before unstable 3-body interactions 
begin. 

The actual enhancement of the gravitational radiation 
rate of course depends on the amount of time spent at high 
eccentricity, but one may place an upper limit on the im- 
portance of Kozai oscillations by computing the radiation 
timescale if the inner binary spends all of its time at e max . 
The orbit-averaged power radiated by gravitational radia- 
tion is given by 
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eccentricity. In Fig. 2 we plot contours of the gravitational 
radiation time r gw — \E / (dE / dt) gw \ in the a-i plane by 



putting e m ax into equation (|10[) . for an equal-mass 6 x 10 
M© binary. This may seem like a gross overestimate of the 
gravitational radiation rate, especially since the shape of 
the Kozai oscillations is in fact such that the binary spends 
more time near e m in than near emax- However since r gw is 
so strongly dominated by periapsis passages at e Ri e max , 
the shift in the co ntours for a realis tic high-e duty cycle is 
only modest. See l|Blaes et al.ll2002T ) for comparison with a 
detailed study of radiation enhancement by Kozai oscilla- 
tions in binaries with initial r gw ~ 10 12 yrs. For a binary at 
a-hard, Kozai oscillations can induce coalescence within 10 10 
yrs in <;20% of cases assuming cosi is uniformly distributed. 
In the remainder of cases the inner binary may survive un- 
til the outer binary shrinks to the point of unstable 3-body 
interactions. 




Figure 2. Upper limit on the importance of Kozai oscillations 
in enhancing gravitational radiation by the inner binary (mi = 
m,2 = 3 X 10 s Mq). The cosine of the initial inclination angle is 
plotted on the horizontal axis, and the inner binary semi-major 
axis is plotted on the vertical axis. Contours are plotted for grav- 
itational radiation timescales of r gw = 10 12 , 10 10 , 10 s , and 10 6 
yrs, if the binary were to stay at maximum eccentricity through- 
out the whole oscillation cycle. The horizontal lines indicate the 
hardening radius and the separation such that a circular binary 
would coalesce on a 10 12 yr timescale. 



1.3 Close 3-body encounters 

If the intruder comes close enough before it causes sufficient 
hardening of the (inner) binary, then a strong 3-body en- 
counter takes place. Strong encounters are characterized by 
a significant transfer of energy between the binary's inter- 
nal degrees of freedom and the COM motion of the binary 
and third body. When the intruder is slow relative to the 
binary's orbital speed vt, in , energy typically flows from the 
inner binary to the outer components, so that the binary is 
more strongly bound after the encounter. This is one mani- 
festation of the negative specific heat characteristic of gravi- 
tationally bound systems. The encounter ends in the escape 
of one of the three bodies, usually the lightest, from the 
system at a speed comparable to Vbin- 

When the lightest body 7713 escapes, momentum conser- 
vation requires that the binary COM recoil in the opposite 
direction with a speed smaller by a factor 7723/ (mi + 7712). 
It is instructive to compare the expected ejection veloci- 
ties of the binary and 7713 with the typical galactic escape 
velocity. For a circular binary with mi = 7712 = mb 4n /2, 
the binding energy at the hardening radius is EB,hard = 
Gml n /8a hard » 6.8 x 10 55 [m Mn /(10 8 Mq)] 3 ' 2 erg. The 
binding energy at the radius where r gw = 10 8 yrs is Es.gw = 
Gml in /8a gw « 6.2 x 10 57 [m bin /(10 8 M Q )] 5/4 erg. The mean 
energy AE harvested from the binary in close encounters 
with slow intruders is about OAEb, though the median AE 
is somewhat lower (|Hills fc Fullertonlll980l) . Energy conser- 
vation implies that the escaper leaves the system with ki- 
netic energy KE s i ng = AE/[1 + 7773/(7771 + 7712)] while the 
binary leaves with KEun.cm = AE/[1 + (mi + m2)/ms] 
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in the system COM frame. For an equal mass binary with 
rribin — 5 x 10 8 Mq, this gives ejection velocities of v S i ng ~ 
290km/s and Vbin ~ 140km/s for the binary at ahard, and 
Wsing ~ 4000km/s and Vbin ~ 2000km/s for the binary at 

Any nonzero eccentricity of the binary will increase the 
semi-major axis corresponding to a fixed r gm , lowering the 
ejection velocities for the binary at a gw . Also any deviation 
from equal masses will result in a smaller fraction of the ex- 
tracted energy being apportioned to the binary and a smaller 
binary recoil velocity. The typical escape velocity for galaxies 
hosting 5 x 10 s Mq BHs is around 1500 km/s, accounting 
for both the stars and the dark matter. From these num- 
bers, it appears that single escapes will be fairly common as 
repeated encounters harden the binary to ~ a gw . However 
accounting for realistic mass ratios and eccentricities (the 
first 3-body encounter tends to thermalize the eccentricity 
even if it starts off circular), binary escapes should be rare. 
Since the binary must come near the escape velocity to re- 
main outside the nucleus for a significant amount of time, 
we do not expect triple interactions to empty many nuclei of 
BHs. We will quantify these statements with our triple-BH 
simulations. 

The formation of triple SMBH systems through inspiral 
of a merging satellite leads to a rather specific initial con- 
figuration. The three BHs start off as a bound "hierarchical 
triple," consisting of an inner binary with a.i n ~ ahard and 
a more widely separated outer binary with semi-major axis 
a.out- For very large a out /ai„ we expect hierarchical triples to 
exhibit very regular behavior; in this case the third body sees 
the inner binary as a point mass and the system essentially 
consists of two independent (inner and outer) binaries. How- 
ever as a ou t/a,i n approaches unity, secular evolution gives 
way to chaotic 3-body interactions in which the orbits di- 
verge and th e system becomes su b ject t o escape of one its 
components. iMardling fc Aarsethl feoOll ) derive a criterion 
for the stability of 3-body systems based on an analogue 
with the problem of binary tides. The most distant intruder 
orbit at which unstable interactions can begin is reliably 
estimated by 
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where R° ut is the periapsis separation of the outer bi- 
nary, a,i n is the semimajor axis of the inner binary, q out = 
W3/(mi +7TI2) is the outer binary mass ratio, and e ou t is its 
eccentricity. This criterion has great practical importance 
due to the high numerical cost of unnecessarily following 
weak hierarchical systems. It specifies an optimal starting 
point for our simulations, which aim to study strong inter- 
actions in 3-body systems starting off as hierarchical triples. 

Naively one might expect a strong 3-body encounter fol- 
lowing a merger with a galaxy hosting a binary, so long as 
the intruder does not induce coalescence of the binary before 
it reaches the center. However the stability criterion implies 
a condition for close interactions much more stringent than 
this. To undergo a chaotic encounter with the inner binary, 
the intruder must reach the stability boundary before the 
outer binary hardens and stalls. A triple system covers some- 
what more stellar phase space than a binary of the same size, 
but not by much for a stable hierarchical system. This means 
that the merger process cannot cause the binary to harden 



by more than around an e-folding for a nearly circular, 
equal-mass system before the intruder arrives at the center. 
Though the order-of-magnitude estimates in the previous 
section make this plausible, further study is needed to deter- 
mine the likelihood of unstable triple interactions in realistic 
merger situations. An eccentric outer binary relaxes the cri- 
terion somewhat, but dynamical friction tends to circularize 
the orbits of satellites with m oderate initial eccentricities be - 
fore they reach the nucleus (|Milosavlievic fc Merrittll200ll ). 
We therefore assume near-circular initial orbits and begin 
each simulation from a weakly hierarchical configuration. 

1.4 Previous work and goals of this study 

Tr iple SMBH systems in galactic nuclei were first considered 
bv lSaslaw et al.l (1 19741 ). who computed an extensive series of 
Newtonian 3- and 4-body orbits, and compared the slingshot 
ejection statis tics to the obser ved structure of extragalactic 
radio sources. IValtonenl l|l976t) included a gravitational radi- 
ation drag force in the 3-body dynamics. He showed that this 
perturbing force could in some cases yield much higher ejec- 
tion velocities than would be possible in Newtonian gravity, 
with associated bursts of gravitational waves. 

The more complex problem of three or four SMBHs 
coming together in the hierarchical merging process and 
interacting in a galactic p oten tial was first address ed by 
iMikkola fc Valtonenl i| 19901 ) and lValtonen et al l (|l994h . who 
experimented wit h a variety of initial BH configurations. 
HcinamakJ |200lf) studied binary-binary scattering in galac- 
tic nuclei using initial conditions (ICs) ba s ed on Extended 
Press S chechter theory (|Lacev fc Coldfl993T ). IVolonteri et al.l 
l|2003af ) followed the formation of triple BH systems in halo 
merger trees tracking the hierarchical buildup of SMBHs 
from ~ 150 Mq seeds in high-cr peaks at z ~ 20. Using a 
simple analytic prescription for the ejection velocities, they 
inferred the presence of a large population of SMBHs and 
IMBHs wande ring through the halo s of galaxies and inter- 
galactic space. Ilwasawa' et all l|2005l ) performed the first full 
N-body simulations of equal-mass triple BH systems embed- 
ded in stellar bulges, an important contribution to our un- 
derstanding of galactic nuclei. Because of the large computa- 
tion time required for each run, they could not statistically 
sample the highly varied outcomes of the 3-body encounters 
as the previous authors did. 

In this paper we study the dynamics of repeated triple- 
SMBH interactions in galactic nuclei. Between close encoun- 
ters we follow the wandering BHs through the galaxy as 
their orbits decay by dynamical friction. We use physically- 
motivated initial BH configurations and mass distribu- 
tions, and updated galactic models characteristic of the 
low-density, massive elliptical galaxies in which SMBH bi- 
naries are most likely stall. We include both a stellar 
and a dark matter component, with the stellar spheroid 
fixed to l ie on the observed m h — a and m B h — M ou i ae 
relations (iTremaine et al.l |2002|; iFerrarese fc Merrittl |2000| ; 



iMarconi fc Huntl2003l ; lMagorrian et alj|l998h . The close en- 
counters are treated using a KS-regula rized Bulirsch-Stoer 
integ r ator p rovided by Sverre Aarseth l|Mikkola fc Aarsethl 
Il990l . ll993l ). The inner density profile is updated through- 
out the simulations to roughly account for core heating by 
dynamical friction and stellar mass ejection. Gravitational 
radiation losses are modelled as a drag force determined by 
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the relative coordinates and velocities of each pair. Each 
simulation takes only a few minutes to run, so we can try 
a variety of distributions of ICs and statistically sample the 
outcomes for each. We use this algorithm to study a vari- 
ety of consequences of the ongoing encounters, such as the 
merging efficiency of BH pairs, the time spent wandering at 
various distances from the galactic center, the distribution of 
final sizes and eccentricities of the binaries remaining in the 
galaxy after a steady state has been reached, and the extent 
of the core scouring caused by the triple SMBH systems. 

Aside from the motivating order-of-magnitude calcula- 
tions in previous sections, this paper does not address the 
question of whether close triple SMBH systems form in galac- 
tic nuclei. We start our simulations from a state that the 
system must reach shortly before the onset of unstable 3- 
body interactions assuming that they occur, and proceed to 
derive the subsequent evolution. Our results may be used to 
argue for or against the occurrence of triple systems in real 
galaxies, as observations support or disfavor the signatures 
that we derive. 

In §2 we describe our model and code methods. In §3 we 
present the results of our study, and in §4 we discuss these 
results and conclude. 



2 MODEL AND METHODS 

2.1 BH mass distribution and halo model 

To get a physically motivated distribution of BH mass ra- 
tios, we associate the formation of the inner and outer bi- 
naries with the last two major mergers in the history of the 
galactic halo hosting t he triple system. W e use Extended 
Press-Schechter theory |Lacev fc C ole 1993) to calculate the 
probability distributions of the halo formation times and 
progenitor masses, and randomly select the parameters of 
the previous two mergers from these distributions. We then 
assign a BH to each progenitor halo using a simple pre- 
scription based on the assumption of a flat galactic rotation 



lLacev fc Coll l| 1993ft derive the instantaneous halo 
merger rate, 



r LC (Mi,M f ,t) 

\da/dM\ Mf exp 



d 2 p 



dM 2 dt 



2 5 C 



7T D{z) 



So 

5c 



— S J— 

' 2D 2 (z) 



a' 2 {M f ) <y' 2 (Mi) 



a 2 {M f ) [l-ff^M^/cr^Mi)] 3 / 2 • ^ 

This equation gives the probability, per unit time per unit 
mass of M2 , of a given halo of mass Mi merging with another 
halo of mass M2 to form a product of mass Mj = Mi + M2 
at time t. Here a 2 (M) is the present-day variance of the 
linear density field on mass scale M, 

a 2 (M) = ——3 / P(k)W 2 {kr)4nk 2 dk, (13) 

where P(k) is the power spectrum of density fluctuations 
today, W is a tophat window function, and r is related to 
M through M = (4/3)7rr 3 p m , the volume times the present- 
day matter density. P(k) is related to the primordial power 
spectrum through the transfer function T{k), which encap- 
sulates the suppression of perturbations on small scales due 
to radiation pressure and damping over the history of the 



u niverse. For T(k) we ado pt the standard fitting formulae 
of lEisenstein fc Hul (| 1998T ) . For the linear growth function 
D(z) we use the approximation 

(5/2) Q m (z)D (z) 



D(z) 



4/7 



(*)-n A (*)+[l + 2*£i] [l + 



Mi 

70 



(14) 



good to within a few perce nt for all plausible values of Q m 
and Oa (|Carroll et al.ll 19921 ). D (z) = 1/(1+2) is the growth 
function for an Einstein-de Sitter universe, Sl m (z) = f2 m (l + 
z) 3 /[f2 m (l + z) 3 + Qa] is the matter density (normalized to 
the critical density) as a function of redshift, and we take 
£Ia(z) = 1 — Q, m (z) assuming the rest of the density is in 
the form of a cosmological constant. S c has the weak redshift 
dependence l|Kitavama fc Sutoll 19961) 

. 3(12rr) 2 / 3 



20 



-[1 + 0.0123 log 10 n m (z)]. 



(15) 



We adopt the cosmological parameters obtained from 
three years of data collection by the Wilkinson Microwave 
Anisotropy Probe (WMAP), Q m h 2 = 0.127, Q, b h 2 = 0.02 23, 
h = 0.73, cr 8 = 0.74, and n s = 0.951 jSpergel et al.ll2006l ). 

Since the merger rate (|12[1 diverges as M2/M1 — > 0, ap- 
plications of the formula that track individual merging halos 
must employ a cutoff mass ratio M2/M1 = A m , such that all 
mergers below A m are treated as smooth accretion rather 
than as discrete mergers (see iManrique fc Salavador-Solel 
1996 for further discussion). The instantaneous rate of ac- 
cretion onto a halo of mass M at redshift z is 

r a (M, t) = / (M' - M)r LC (M, M' , t)dM' . (16) 



To get the growth history ("accretion track") of a halo 
of mass Mo at time to due to accretion since the last 
merger, one need only solve the differential equation dM/dt 
= r a [M(t),t], subject to the initial condition M(to) = Mo. 
We integrate this equation backward in time using a 4th- 
order Runge-Kutta method to get the accretion tracks of the 
halos in our simulations. Since we are interested in BH bi- 
nary formation, we loosely associate A m with the halo mass 
ratio such that tidal stripping of the satellite would prevent 
the eventual merging of the two nuclei. N-body simulations 
of galaxy mergers place this mass ratio in the range A m ~ 
0.1 — 0.3, depe nding on the density and orbital par ameters 
of the satellite (iTaffoni et al.ll2003l ; lcolpi et al.ll 19991 ) . Hence 
our canonical choice is A m = 0.3, and we also try values of 
A m = 0.1 (runs Dl) and 0.5 (runs D5), the latter being the 
halo mass that corresponds to a stellar mass ratio of ~ 3:1 

in our prescrip tion. 

Following ISalvador-Sole. Solanes. fc Manriquel l|l998l) . 
we write the probability, per unit time, of a halo with mass 
Mf at time t arising from a merger with a smaller halo of 
mass between M and M + dM (the "capture rate" ) as 

r c (M,M f ,t) dM = 

r LC (M,M f ,t)e[M f - M(l + A m )] ^ff'*]) dM ' ( 17 ) 

the EPS merger rate excluding halos below the threshold 
A m , and weighted by the number of mass M halos per unit 
halo of mass Mf. 

The rate at which halos of mass M/ form through all 
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mergers at time t is 

1 r M f 

f/(M/,*)w- / r c (M,M f ,t)dM. (18) 

z J M f A m /(l+A m ) 

The probability distribution function (PDF) of formation 
times of halos with mass Mo at time to is 

$/(M ,t) =r f [M(t),t]e- f *° r f lMit " ) ' t ' ]dt '. (19) 
Given a formation time t/ and corresponding mass M(t/) 
along the past accretion track of Mo, the mass of the larger 
progenitor Mi is distributed according to 

2G(Mi,M) 



*„[M(t/),Afi 



r M/(l + A m ) 



+A„ 
/(1+A„ 



G(M',M)dM' 



where 



G(M',M) 



\da(M')/dM' 



\M) 



1 -3/2 



2 (M') 



(20) 



(21) 



M'(7 2 (M') 

By choosing formation times and progenitor masses ran- 
domly according to (I19|) and (|20p . we capture the 
stochasticity of the intervals between mergers above 
A m , but treat merging below this threshold only in 



the mean. See Raig, Gonzalcz-Casado, & Salvado r-Sold 



ll200ll). ISalvador-Sole, Solanes, fc Manriqud (|l998l ). and 
iManrique fe Salavador-Sold ( 19961 ) for further details and 



derivations of (|19p and (|2U|) . Fig. 3 shows the distribution of 
formation times, progenitor masses, and accretion tracks for 
a present-day 5 x 10 13 M Q halo for A m = 0.1, 0.3, and 0.5, 
and the accretion tracks for 1, 2.3, and 5 xlO 13 Mq halos 
with A m fixed at 0.3. All accretion tracks are normalized 
to the present-day mass Mo. Note the insensitivity of the 
shape of these tracks to Mo, as expected for masses above 
the critical mass M* ?s2x 1O 12 M . 

Our algorithm for generating the BH masses is illus- 
trated schematically in Fig. 4. For each run we begin with a 
halo of mass Mo = 5 x 10 13 Mq at time to = t(z — 0), choose 
it's formation time tfo randomly according to equation (119[) . 
and find the mass M/o = M(tfo) along its accretion track 
at that time. The mass M/o is assigned to the dark matter 
halo hosting the triple BH system, and the physical time for 
the run to end if other termination conditions are not met 
first is set to to — tfo- To explore the dependence of the re- 
sults on the absolute mass scale, we also try beginning with 
a 1 x 10 13 Mq halo (runs HI). 

We model the halo as a Hernquist profile 
jHernquistj Il990h. which is identical to an NFW pro- 
file I Navarro et al.| [l997) in its inner regions if the scale 
radius an is related to the NFW scale radius by an = 
OiVFW\/2[log(l + c) — c/(l + c)], where c is the halo con- 
centration defined by ajvFW = r V i r /c. The H ernquist model 



falls off as r 4 instead of r 3 
120051 ). The virial radius r v 



far outside an ( Spr ingel et al.l 
(M, z) is given by 



r vir (M, z) = 
364 



1 + z 



M ha l fl m (z)l8TT 2 



21 1/3 



w i3 h- i MQ n„ 



h~\pc, (22) 



where A c 



18t t 2 + 82\Q m (z) - 1] - 39[fi m (z) - l] 2 



i|Barkana fc Loebl 120011 ). and c roughly fol lows the me- 
dian r elation from the ACDM simulations of iBullock et al.l 
(l200ll) . c w 9.0[(2.1 x 10 13 M Q )/M ha , o ] 013 /(l + z). The z 
dependence of r„i r and c nearly cancel to make oh depend 
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Figure 3. Upper left: Probability distribution function (PDF) of 
formation redshifts given by equation (1 1 Q I t for a 5 X 1O 13 M0 halo 
at 2 = 0. Upper right: PDF of masses of the larger progenitor of 
the same halo given by equation H20I I. normalized to the mass of 
the merger product, Mf = 5 X 1O 13 M0. Lower left: Past accretion 
tracks of a present-day 5 X 10 13 Mq halo back to z = 3, normalized 
to the mass at z = 0. Lower right: Normalized accretion tracks 
for three different halo masses. 



only weakly on redshift, so we simply use the z = relation 
between M^aio and an in our simulations. 

We choose the mass Mi of the larger progenitor of M/o 
randomly according to equation (|20|l , and assign a mass M2 
= M/o — Mi to the smaller progenitor. Before the merger 
the larger progenitor is assumed to have hosted a BH binary, 
while the smaller one hosted a single BH. Repeating the 
procedure used for Mo, we assign formation times t/i and 
t/2 to Mi and M2 using equation (|19[) . and choose progenitor 
masses Ma, M12, M21, and A/22 according to equation <|20l) . 

Having constructed a set of progenitor halos, we now 
need a BH-halo relation mbh{Mhalo, z) to complete our al- 
gorithm. We obtain such a relation by equating the halo 
virial velocity v V i T to the circular velocity v c of the stellar 
spheroid, and using empirical v c — a and a — nibh corre- 
lations to connect Vc to mn,, similar to t h e app roaches in 
Eric kcek et al, I (|2006l ) and lWvithe fc Loebl |2005). Combin- 
ing 



; 343 VI + z x 

/ M ha l 

VlO^/t-iMe 



1/3 



£ 2.1) 



A c 



1/6 



n m (z)i8^\ kms " (23) 

|Barkana fc Loebl l200lf) with v c » 314[a/(208 km/s)] 84 
km/s dFerraresd 120021) and o7f208 km/a) w [m bh /(1.56 x 
10 8 M^)] 1/4 02 (jTremaine et al.ll2002l ). we arrive at the rela- 
tion 



M h 



1O 12 M 



= 8.28 



M b 



(lO 8 W. 



7(2) 



(24) 



where j(z) = (1 + z)- :1 / 2 [(n m /n m (z))(A c /l87v 2 )}-^ 2 . 

In our canonical runs we set the masses of the in- 
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Figure 4. Schematic diagram of our algorithm for generating 
the BH mass distribution. First we select a formation time zjq 
for the halo Mo hosting the BH triple system randomly from 
equation 1191 . Given Zfo and M(Zfo), we select two progenitor 
masses Mi and M2 according to equation ||20|I . assign the binary 
to the larger one and the third BH to the smaller one. We repeat 
this process going back one step further in the "merger tree" to 
get the masses of the binary constituents. 



ner binary members to mbh(Mu, z/i) and mbh(Mi2, Zfi), 
and that of the intruding BH to m^^M^i + -^22,2/2)- 
Note that in this prescription the intruder is usually lighter 
than the heavier binary member, so that most of the 3- 
body interactions result in an exchange. To examine the 
effect of more interactions without exchange, we try choos- 
ing ?7if,fe[max.(Af2i; M22), Zfi] for the intruder mass in runs 
(MX). As there is neither a direct causal relationship be- 
tween rribh and Mhaio predicted by theory (|Wvithe fc Loebl 
120051 ) nor a tight correlation directly observed between these 
two variables, and we know that identical halos may host 
galaxies of different morphologies and occupation numbers, 
mbh{Mhaio, z) should be taken with something of a grain of 
salt. Nevertheless it is a useful way to generate simple but 
physically-motivated BH mass distibutions when no infor- 
mation other than the halo mass is available. 

We make one final modification to the set of BH masses 
used in our simulations. If the outer binary's hardening 
radius lies outside the stability boundary given by equa- 
tion (|f 1 1 with din = cihard, then the decay of the outer orbit 
is expected to stall before a strong encounter can begin. To 
roughly account for this we exclude all ICs where f/, a ut > 
3/iin- The final distribution of BH mass ratios is shown in 
Fig. 5 for A m =0.1, 0.3, and 0.5. In the upper panel we plot 
the inner binary mass ratios, while the lower panel shows 
the distribution of mfc.i n /m eS c, where m e3C is the mass of 
the lightest BH and mun is the sum of the masses of the 
other two BHs. This ratio determines the binary recoil speed 
when the lightest BH is ejected from the system. The total 
BH mass is typically ~ 6 x 10 8 Mq in our canonical runs. 
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Figure 5. Distribution of BH binary mass ratios. Upper panel: 
Inner binary mass ratio mi/m2 > 1. Lower panel: m bin /m esC; 
where m eS c is the mass of the lightest BH and m|,i„ is the sum 
of the masses of the other two BHs. 



2.2 Stellar spheroid model 

To complete the galactic model we surround the BH sys- 
tem by a smooth stellar potential superimposed on the dark 
m atter halo. The s tars a re modelled using the "^-models" 
of iTremaine et ail {l994), with a sharp break to shallower 
slope —7 added at rb <C a: 

77 Ma _ 

!>{l) = < ^r*-v(r + a)^- Pr > {r) ' * r > n > (25 ) 

if r < r b . 



47T r 3 - 7 !^ + a) 1+r > 
i Pr ? (r b )(r/r !) )- 7 

Our canonical model is the r\ = 2 (Hernquist) profile, and we 
also try r\ = 1.5 (runs SC) to explore the effect of a steeper 
inner profile and higher central density (p ~ 800 Mq pc -3 
for -q — 1.5 vs. p ~ 180 Mq pc -3 for the Hernquist profile 
at the BH radius of influence). Tb and 7 were initialized 
to reflect the cusp destruction caused by the inspiralling 
BHs in reaching their initial configuration, and were updated 
throughout the simulation to account for the continued core 
heating and mass ejection. Our algorithm for updating the 
core is described further in §3.5. 

The parameters M and a in the ry-models were set 
based on the tight correlations ob served between SMBH 
mass mbh and stellar bulge mass lIMagorrian et al. I Il998l ; 
iMarconi &; Huntll2003l; IPeng et al.l [20051) and velocity dis- 
persion dFerrarese fc MerrittJ |200d; iGebhardt et al.1 bood ; 
ITremaine et all I2002T ). IMarconi fc Hunt! (|2003l ) found the 
relation M M ho3 = (4.06 x lO lo M0)(m^/lO 8 M ) 104 be- 
tween rribh and the virial mass M V i r — kR e a^/G of the 
stellar bulge, where R e is the half-light radius and <r e is 
the effective bulge velocity dispersion. They set k — 3 (k 
would be 8/3 for an isothermal sphere) to get an average 
ratio of unity between M V i r and the dynamically measured 
masses Md yn of galax ies with more direct s tellar-dynamical 
mass determinations l|Gebhardt et al.| [2003'l. a e is typically 
measured over either a c ircular aperture of radius R e /8 
dFerrarese fc Merritulioooh or a linear aperture out to i? e 
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jCebhardt et alj|2000l ) - iTremaine et all (|2002h discuss the 
essential agreement between the velocity dispersions mea- 
sured in these two ways. Thus for each model (77 = 2, 1.5) 
we compute the projected radius R e = Kia containing half 
the integrated surface brightness (assuming a constant mass- 
to- light ratio), and velocity dispersion <j\ = K 2 GM/a at ra- 
dius R e / 8. Here ki and k 2 are constant s depending on the 
density profile (see ITremaine et al.ll 19941 for relevant formu- 
lae). The parameter M is then chosen to satisfy 3R e a^/G = 
3kiK 2 M = M M H03,{m bh ) M = A/ A /iW(3«:iK2). For the 
Hernquist model with «i = 1.815 and k 2 = 0.104, Mh = 
1.76M A fH03 = (7.15 x lO lo M )(m6h/lO 8 M o ) 104 . The scale 
radius is then obtained from a — GMMH03/3Kial h (m,bh), 
where abh(rnbh) is the velocity dispersi o n com puted from 
the rribh — a relation of ITremaine et all (|2002T h In each of 
these relations rribh is set to the total mass of the triple BH 
system. 

In a perfectly smooth, spherically symmetric galactic 
potential, BHs ejected on distant radial orbits return di- 
rectly to the center to interact strongly with any other nu- 
clear black holes. Since real galaxies are clumpy and triax- 
ial, the interaction will more realistically be delayed until 
the orbit of the ejected BH decays by dynamical friction. To 
mitigate this problem we flattened the 77-models by adding 
two low -order spherical harmonic terms to the spherical po- 
tential (|de Zeeuw fc Carolldfl996h : 

V{r) = u(r) - v{r)Y 2 °(9) + w(r)Y 2 2 (6, 0), (26) 

where u(r) is the potential of the spherical 77- model, v(r) = 
-GMnr*!- 1 /( r + r 2 ) v+1 , and w(r) = —GMr' i r v ' 1 /(r + 
ri) rl+1 . Since near-sphericity in the inner regions is prob- 
ably a necessary prerequisite for the survival of the inner 
binary for of order a Hubble time until the next merger 
|Merritt fc Poonll2004l ; iBerczik et al.ll2006l ). the parameters 
n,...,r4 were chosen to give a spherical profile near the 
galactic center, and axis ratios approaching 1.3 and 1.5 for 
r 3> a. A similar triaxial modification was applied to the 
dark matter halo, and the relative orientation of the stel- 
lar and halo potentials was chosen randomly. By misalign- 
ing their axes we eliminate any artificial stable orbits (e.g. 
along the long axis of an ellipsoid) near which ejected BHs 
tend to return on a perfectly radial orbit to the center. This 
triaxial modification had the desired effect of preventing fre- 
quent strong encounters at periapsis on distant orbits, but 
had little influence on the global outcome statistics. 



2.3 Initial BH configuration 

We assume that the 3-BH system starts off as a hierarchical 
triple on the verge of unstable 3-body interactions. In our 
canonical runs we initialize the inner binary semi-major axis 
a,i n to ahard- To study the effect of varying a,i n we also try 
runs with a in — 3rh (runs BA) and a in = rh/3 (runs SA). 
The outer binary semi-major axi s a ou t is set by the stability 
criterion of iMardling fc Aarsethl (|200ll ). equation fTT]). The 
initial eccentricity of the inner (outer) binary was chosen 
uniformly between 0.0 and 0.2 (0.3), in accordance with the 
low eccentricities found in galaxy merger simulations where 
dynamical friction tends to circularize the orbits of sate l- 
lites as they spiral inward (|Milosavlievic fc Merrittll200ll ). 
The three Euler angles of the intruder's orbital plane were 
chosen randomly relative to the reference plane of the binary 



orbit, as was the phase of the initial periapsis of the binary. 
Both orbits were always started at periapsis; since many or- 
bital periods elapse before unstable interactions begin, the 
relative phase is effectively randomized in any case. Having 
defined an initial configuration of three BHs embedded in 
a stellar+dark matter potential, we next describe how we 
evolve the system forward in time. 



2.4 Code method 

We treat the close 3-body encounters using Sverre Aarseth's 
Chain code, an implementation of the N-body regulariza- 
tion t echniq ue of Mikkola and Aarseth l|Mikkola fc Aarsethl 
Il99d . Il993t ). The masses are first ordered so that neigh- 
bors in the chain are the do minant two-body interactions , 
then the KS-transformation (|Kustaanheimo fc StiefelHl965T ) 
is applied to neighboring pairs. This transformation elim- 
inates the singularity at r — > in Newtonian gravity and 
transforms the equations of Kepl erian motion to the sim - 
ple harmonic oscillator equation (Sticfcl & Schcifclc 1971). 
External perturbing forces of arbitrary strength depending 
on the coordinates, velocities, and/or time are simply incor- 
porated into the formulation (though of course singularities 
in these perturbing forces need not be eliminated by the 
change of variables). We use this to add a galactic poten- 
tial (§2.1-2), a gravitational radiation back-reaction force, 
and a stellar-dynamical friction force on the intruding BH. 
The regularized equations of motion are integrated us ing the 
Bulirsch-Stoer (BS) method jBulirsch fc Stoerlll966h based 
on Romberg extrapolation. For unperturbed sinusoidal mo- 
tion, the BS integrator requires only two or three timesteps 
per orbital period! 

When the binary and third body are far apart we switch 
to two-body motion (of the single BH and binary COM) us- 
ing a 4th-order Runge-Kutta (RK4) method. We simultane- 
ously evolve the binary semi-major axis and eccentricity us- 
ing orbit averaged equations, da = [(da/dt) st + (da/dt) gw ]dt 
and de = [(de/dt) gm ]dt, where (d/dt) at and (d/dt) gw are 
the contributions from stellar interactions and gravitational 
radiation. The timesteps are adaptively controlled with a 
simple step-doubling scheme: at each step the 14 numbers 
{xi, xg; hi, v§\ a, e} are all required to remain the same 
to within an error e = 10 _n under doubling of the step size. 
To avoid wasting computation time when any of these values 
approach zero, we accept agreement to n decimal places as 
an alternative criterion for convergence. For the calculations 
reported in this paper we set n = 12. 

The relative perturbation to the binary from the third 
body at apoapsis, 

5F= .f^V (27) 
mm(mi, 7712)0 

is used to decide which integration method to use at any 
given time. Here r ap is the apoapsis distance between mi 
and m2, 7773 is the intruder mass, and d is the distance of 
the intruder from the binary COM. We switch to two-body 
RK4 integration each time SF falls below 5 x 10 -5 and call 
the Chain code again when 8F reaches 5 x 10 -4 . We choose 
different 8F thresholds for beginning and ending close en- 
counters to prevent overly frequent toggling between the two 
methods. When < 3 BHs remain in the simulation (after co- 
alescence of the inner binary or escape of one or more BHs 
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from the galaxy), we primarily use the RK4 integrator, but 
call the Chain code to treat very close two-body encounters. 
Since chain regularization is defined only for three or more 
bodies, we add a light and distant "dummy" particle when 
using this method for two-body motion. 

During the two-body motion we declare the single BH 
or binary (remnant) to have escaped if its distance from 
the galactic center exceeds 500 kpc and its specific energy 
E = $(r, 9, <f>) + iv 2 exceeds E esc , the energy needed to es- 
cape from r — to infinity. We declare the binary to have co- 
alesced during a close encounter when (i) ri2 < 3r s t, where 
r s t is the Schwarzchild radius of the larger member of the 
pair, or (ii) \a/a\ gr < 0.1td yn and \a/a\ gr < 50 yrs while 
SF < 10~ 3 , where tdyn is the current outer binary dynami- 
cal time. During the RK4 integration we require that \a/a\ gr 
< 50 yrs or r\i < l.l(r s i +r S 2) at periapsis, where r s i,2 are 
the Schwarzchild radii of the two binary members. Upon 
coalescence we replace the pair with a single body of mass 
mun an d the COM position and velocity. A run ends when 
(a) only one SMBH remains in the galaxy and it has settled 
to the center of the potential by dynamical friction; (b) two 
BHs remain and have formed a hard binary at the galactic 



center; (c) the physical time exceeds t„ 



= to -t 



/ti- 



the 



current age of the universe minus the halo formation time; 
(d) all BHs have escaped the galaxy; or (e) the physical time 
spent in a call to Chain exceeds a maximum allowed time 
t c hn ■ The last condition is added to avoid spending too much 
computation time on very long close encounters. 



2.4-1 Treatment of stellar- dynamical friction 

During the two-body evolution we apply a dynam- 
ical friction force given by Chandrasekhar's formula 



(28) 



jChandrasekharll 19431 ; binney fc Tremaindfl987h . 

'<Fv\ _ 4ivG 2 pm In A [erf (X) - 2Xe~ x2 /0F] . 
dt) dJ v 2 

where X = v/(y2a), to the single BH and binary COM. 
The factor in square brackets ~ 1 for v 3> o and ~ 0.75X 3 
for v <C o. We take 

In A = max (in , i\ (29) 

[ [ Gm J J 

for the Coulomb logarithm, where r is the BH's distance 
from the galactic center. For p in equation (|28p we use 
min[p(r), p(ri„f)], effectively capping the density at its value 
at the BH radius of influence, rt„f — Gm/a 2 , when the BHs 
pass through the core. 

The semi-major axis a of the binary also evolves under 
stellar-dynamical friction as it wanders through the galaxy. 
However Chandrasekhar's formula applied separately to the 
binary constituents does not give a good description of this 
evolution, since the hard binary loses energy through close 
3-body encounters with stars, while equation (|28p relies on 
the assumption that the energy loss is dominated by weak 
two-body encounters. We approximate the evolution of a 
using a formulation for the decay rate of a hard, massive 
binary in a uniform a nd iso tropi c sea of s t ars de veloped in 
iMikkola fc Valtonenl |l992l) and lOuinlanl (|l996l ). The for- 
mulation was calibrated with an extens ive series of 3-body 
scattering experiments inlQuinlanl (Il996 | ) and teste d against 
N-body simulations in lMikkola fc Valtonenl (|l992T l. The bi- 



nary decay rate is given by 
da s 
IE 



pHa 2 



(30) 



where the hardening rate H can be ap proximated by the 
empirical fitting function (|Quinlanlll99q ) 

(3D 



H ; 



[1 + [o/wf] 



41 1/2 ' 



Here w — 0.85^/G min(mi,m2)/a is the characteristic ve- 
locity distinguishing the hard binary regime - stars with 
v J>™ cannot be easily captured into bound orbits and pref- 
erentially harden the binary in close encounters. In our sim- 
ulations the binary COM is often speeding through the stel- 
lar medium at v cm J>cr after an energetic ejection, so the 
stellar medium looks "hotter" in its frame of reference. To 
account for this we replace a in equations (|30|l and (|31|) 



with <7* 



V v ' 2 m + °" 2 5 a good approximation since H is 



not very sensitive to the shape of the distribution func- 
tion (e.g. H « 16 for a Maxwellian vs. H ~ 18 for a uni- 
form velocity distribution). For p in equation (|30|l we took 
min[p(r), p(ri n f)] as we did for the drag on the COM. We 
ignored t he mild eccentr icity evolution (de/dt) at , which is 
shown in lOuinlanl (|l996l ) to be far weaker than that pre- 
dicted by Chandrasekhar's formula for hard eccentric bina- 
ries. 

When the amplitude of oscillation of one of the two 
masses falls below Gm/2a 2 , we stop integrating its motion 
and place it at rest at the galactic center until the second 
body returns to within a distance of twice the break radius, 
2rv If the settled mass is the binary, then we also stop up- 
dating its semi-major axis for stellar hardening, assuming 
that it clears out its loss cone and stalls once it stops mov- 
ing about the nucleus and encountering new stars. Since the 
total mass in loss cone stars is small compared to the BH 
mass in the low-density galaxies that we consider, to good 
approximation the binary stalls as soon as the replenishing 
mechanism (motion) shuts off. 

During close encounters between the three BHs an 
orbit-averaged prescription for stellar-dynamical friction is 
not feasible. However the triples are still marginally sta- 
ble at the boundary given by equation (JTTJ , so we apply a 
drag force given by Chandrasekhar's formula with In A = 1 
to the intruder at the beginning of each run. At the on- 
set of chaotic interactions in the first encounter (defined 
loosely by the first time the closest pair is not formed by 
the original binary members) this perturbation is shut off, 
and it remains off in all later close encounters. Fortunately 
the chaotic interactions occur on timescales very short com- 
pared to a dynamical friction time, so it is valid to neglect 
stellar dissipation during close encounters. 



2.4-2 Treatment of gravitational radiation 

Gravitational radiation is modelled using the 0[(v/c) 5 ] post- 
Newtonian (2.5PN) back-r eaction acceleration computed by 
iDamour fc Deruellel (|l98lf ). evaluated in the two-body COM 
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frame (e.g. Gultekin et al. 2006), 



dv\ 4G m 2 
dt 5c 5 mi + m 2 



mira 2 



|f(r • v) ^ 



34G(mi + ra 2 ) 



3r 



+ 6u + 



"1- 



6G(rm 



m 2 ) 



- 2v 



']}■ (32) 



r = fi — r 2 and t? = ui — ?7 2 are the relative positions and 
velocities of the two masses. We sum the force linearly over 
all pairs, a valid approximation provided the perturbations 
from the third body and other external tidal forces are in- 
stantaneously small at periapsis. When averag ed over a com - 
plete orbit, equation (|32[) is equivalent to the lPetersl (|l964T l 
equations for the binary semi-major axis and eccentricity, 



da 
~dt 

de 
dt 



64 G 3 mim 2 (mi + m 2 ) 1 + 5j e + § e 



5 c 5 a 3 
304 G 3 mim 2 (mi + m 2 ) 



(l_ e 2)7/2 

e + — e 3 



15 



(!_ e 2)5/2 



(33a) 



(33b) 



However when |t)-fj comes close to one on hyperbolic orbits, 
so that (f ■ v) 2 — > ii 2 , E — F\ ■ v\ + F 2 ■ « 2 = miai ■ « as 
given by equation (I32|l becomes positive, though we know 
physically that gravitational waves can only carry energy 
away from the system. To give the correct answer averaged 
over an orbit, this positive contribution must be cancelled 
by extra energy loss near periapsis, making the equation 
potentially sensitive t o numerical error. This effe ct is much 
less pronounced in the lDamour fc Deruellel (Il98ll ) form than 
in other expressions derived for the radiation back-reaction 
acceleration - they derived the formula specifically for prac- 
ti cal use on the problem of two point masses (see Appendix 
of lLee| [l993 and references therein). 

For computational ease we neglect the lower-order 
1-2PN terms (precession of the periapsis) in the post- 
Newtonian expansion. Though much larger in magnitude 
than the radiation reaction force, these terms are unimpor- 
tant in the statistical sense because they cons erve the in- 
trinsic properties of the sy stem, such as energy l|Kupi et al.l 
2006; I wasawa et al.ll2005h . We need not concern ourselves 
with relativistic precession destroying the Kozai resonance 
since the semi-major axis ratio given by equation (jTTJ is 
much smaller than that of equation ©. 



2.5 Code tests and energy errors 

One way to establish the reliability of our integration meth- 
ods is to test them on problems with known solutions. Fig. 
6 shows an example on the two-body problem with gravi- 
tational radiation. The upper panels show the evolution of 
the semi-major axis o and eccentricity e of four decaying el- 
liptical orbits, computed using (a) our RK4 integrator and 
equation (|32p . with an error tolerance of e = 10 -9 , (b) the 
Chain code a nd equation (|32[) . with e = 10~ 14 , and (c) the 
iPetersI l|l964h equations (133 ^ - In each case the initial semi- 
major axis ao was chosen to give a gravitational radiation 
timescale of \a/a\ « 10 7 yrs, and the four curves (from bot- 
tom to top) are for eccentricities of 0.0, 0.5, 0.9, and 0.99. 
The agreement of the three computation methods demon- 
strates the reliability of both the RK4 integrator and our 
implementation of the Chain code in handling dissipative 
forces. 

The lower panels show two hyperbolic orbits with peri- 
apsis distances around 30 times the Schwarzchild radius r s b 



t (10 7 yrs) 
2 4 6 




Figure 6. Code tests on the Newtonian two-body problem with 
gravitational radiation. Upper panels: Evolution of elliptical or- 
bits computed using RK4 integration of the Peters equations 
(black solid), RK4 integration with the Damour & Deruelle (DD) 
radiation back-reaction acceleration (red dashed), and Chain in- 
tegration with the DD acceleration (diamonds). The ICs were 
chosen to give \a/d\ = 10 7 yrs at the beginning of each integra- 
tion. Left: Semi-major axis evolution for initial eccentricities of eo 
= 0.0, 0.5, 0.9, and 0.99 (bottom to top). Right: Eccentricity evo- 
lution for eo = 0.5, 0.9, and 0.99. The curves' indistinguishability 
demonstrates the reliability of all three methods. Lower panels: 
Hyperbolic orbits with impact parameters b set to 80% and 120% 
of the critical value for gravitational radiation capture, computed 
using the DD acceleration in the Chain code. Left: 0.8b cr i t ; BH is 
captured. Right: 1.2b cr ; t ; BH is not captured. The blue asterisks 
are points along the Newtonian trajectory (without gravitational 
radiation). The deviation from the Newtonian trajectory after pe- 
riapsis can be seen in both plots, even though the energy remains 
positive in the latter. 



of the larger BH, computed using equation (|32|l in Chain. 
The RK4 integrator was found to fail some tests on very 
close approaches from hyperbolic orbits with gravitational 
radiation, so we treat all such approaches using the regular- 
ized Chain code in our runs, even during the unperturbed 
binary evolution. The blue asterisks are points along the 
Newtonian orbits while the red solid lines show the trajec- 
tories with gravitational radiation. There is a simple analytic 
expression for the maximum periapsis distance for gravita- 
tional radiation capture from a hyperbolic orbit, 



85^2ttG 7/2 rmm 2 (mi + m 2 ) 



3/2 



2/7 



iw. ■ (34) 

where mi and m 2 are the masses of the two bodies and Doo 
is their relative velocity at infinity. The orbit on the lower 
left begins at 80% of the critical impact parameter and the 
incoming BH is captured. On the right the intruding BH 
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Figure 7. Energy errors. Upper panel: Errors given by equa- 
tion (1351 . for each code stage: close triple encounters, integrated 
with Chain (red); unperturbed binary evolution with the RK4 
integrator (blue); and very close two-body encounters computed 
with Chain after the binary has coalesced (green). The black 
(heavy) histogram shows the energy errors for close triple en- 
counters normalized to the initial energy instead of the energy 
dissipated. Lower panel: Effective energy errors for the entire run, 
computed from equation ((36} • 



starts at 120% of the critical impact parameter and is not 
captured, though the deviation from the Newtonian trajec- 
tory due to the energy radiated at periapsis can be seen 
on the way out. We tried iterating over impact parameters 
close to the critical value and found that the code reproduces 



equation (|34|) to within a part in 10 for periapsis distances 



30r s , and to within a part in 10 for r pe 



3r s 



We also evaluated the performance of the code by re- 
peating our canonical set of 1005 runs with a static inner 
profile to check the precision of energy conservation. In Fig. 
7 we histogram the energy errors, computed as 



E + 



Vi + Fg 



Vi)dt\ - E f 



Eo - E f 



(35) 



where Eo and Ef are the initial and final energies, and the 
two terms in the sum under the integral are the work done 
by dynamical friction and gravitational radiation during the 
current stage of the code. In the upper panel we separately 
plot the errors for close 3-body encounters, RK4 integra- 
tion of the unperturbed binary motion ("far"), and close 
two-body encounters computed with Chain during the un- 
perturbed binary evolution. The plot includes all code stages 
where the energy dissipated was at least 10 -3 in code units, 
or about a part in 10 5-6 of the initial binding energy of the 
system. The black (heavy) histogram shows the errors for 
close encounters normalized to the initial energy instead of 
the dissipated energy in the denominator of equation (|35p . 
since the energy dissipated was very small in many close en- 
counters. In the lower panel we combine the energy errors 



from the various code stages to get an effective energy error 
for each entire run, 



trun - AE 1+ AE 2 + ... + AE n ■ W 
We had to combine the separate errors to obtain e run since 
the galactic potential is handled slightly differently during 
different stages of the code, e.g. the triaxial modification is 
applied only during the RK4 integration. In a large major- 
ity of cases e run falls between 10 -12 and 10 -9 , and energy 
is conserved to better than a part in 10 4 in every run. The 
excellent energy conservation gives us confidence in the ro- 
bustness of our integration methods. 



3 RESULTS 

3.1 Outcome statistics 

We begin with an overview of the outcomes of our 3-body 
simulations. In subsequent sections we focus on various ef- 
fects in more detail. Our data consists of eight sets of 1005 
runs, each sampling a different distribution of the ICs. A set 
of 1005 runs took anywhere from ~4 to ~30 hours to finish 
on five 2.0 GHz Opteron processors, depending on the ICs. 

In our canonical runs (CN), we chose A m = 0.3 for 
the threshold merger mass ratio, modelled the stellar bulge 
as a Hernquist (rj = 2) profile, started off the inner binary 
at tthard, and generated the ICs from a 5 x 10 13 M© halo 
at z = 0. In each of the remaining runs we varied one of 
these assumptions. Runs Dl and D5 used A m = 0.1 and 
A m = 0.5 to explore the effects of widening or narrowing 
the range of BH mass ratios. In runs MX we assigned a mass 
m6h[max(M2i, M22), z/a] instead of mbh(M 2 i + M22, 2/2) to 
the intruding BH, as discussed in §2.1. In runs BA and SA we 
started off the inner binary at 3ahard and ahard/3 instead 
of at a^ard- We initialized the stellar bulge to an r\ — 1.5 
profile in runs SC, to explore the effect of a steeper inner 
cusp. Finally in runs HI we generated the ICs from a 1 x 10 13 
Mq halo at z = 0, for total BH masses of ~ 5 x 10 7 M , 
about an order of magnitude lower than in our canonical 
runs. Table 1 summarizes the outcomes. 

The first two rows give the percentage of cases in which 
(i) one BH pair coalescenced by the end of the run (i.e. by 
the time since the last major merger) , and (ii) the remaining 
two BHs also coalesced within the run time. At least one pair 
coalesced in a large majority of the runs for each set of ICs 
that we tried. The new system formed from the third BH 
and binary remnant also coalesced in ~10-20% of the cases. 
Since we assume that stellar hardening of the new binary 
shuts off at ahard, it can only coalesce by gravitational radi- 
ation from a highly eccentric orbit; we will discuss this topic 
further in §3.2 and §3.4. The coalescence rate is somewhat 
lower (~ 68%) in set Dl, since (a) the hardening effect of 
the third body is lessened for more extreme mass ratios, and 
(b) mergers with mass ratios as low as A m =0.1 are more 
frequent, so the run time is typically shorter. Naturally the 
coalescence rate is somewhat higher (95%) in runs SA, where 
we begin with a tighter binary (ao = a^ard/3, r gw ~ 10 11-12 
yrs). Coalescence is also significantly less common in runs 
HI. This can be understood in light of equation © in §1.1. 
The separation between the scale set by the stellar kinemat- 
ics (ahard) and that set by gravitational radiation (a gw ) is 
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Table 1. Summary of outcomes for eight different distributions 
of the ICs. See text for explanation of the entries. 



Outcome 


CN 


Dl 


D5 


MX 


BA 


SA 


SC 


HI 


Coalescence 


One pair 


87% 


68% 


89% 


84% 


84% 


95% 


84% 


75% 


Two pairs 


15% 


13% 


18% 


16% 


22% 


13% 


16% 


10% 


Escape 


Single 


15% 


18% 


17% 


19% 


21% 


14% 


14% 


22% 


Double 


0.0% 


0.0% 


0.1% 


0.0% 


0.1% 


0.2% 


0.1% 


0.3% 


Wander 


37% 


51% 


38% 


42% 


42% 


26% 


45% 


50% 


Final state 


Binary 


60% 


67% 


54% 


57% 


51% 


65% 


56% 


63% 


Single 


38% 


31% 


45% 


41% 


47% 


33% 


43% 


35% 


No BH 


1.3% 


2.0% 


1.2% 


1.5% 


1.4% 


1.4% 


1.2% 


1.8% 


Core 


(M def ) 


1.41 


1.20 


1.50 


1.38 


1.31 


1.96 


1.42 


1.61 


AM de/ 


0.48 


0.37 


0.57 


0.44 


0.52 


0.44 


0.57 


0.61 


Termination 


Sing escape 


14% 


16% 


16% 


18% 


18% 


14% 


12% 


21% 


New binary 


61% 


44% 


61% 


55% 


53% 


73% 


54% 


49% 


T.O. far 


22% 


30% 


22% 


23% 


20% 


12% 


31% 


28% 


T.O. Chain 


2.1% 


10% 


1.7% 


4.4% 


8.4% 


0.3% 


1.6% 


1.5% 


Crashed 


0.5% 


0.1% 


0.3% 


0.3% 


0.7% 


0.2% 


0.6% 


0.1% 



proportional to m^J/ 4 . Hence in lower mass systems, coales- 
cence is less likely relative to escape. This observation mo- 
tivates future study of triple BH dynamics in much lighter 
systems. 

The next three rows of the table give the fraction of 
runs in which (i) the single BH escaped the stellar bulge + 
halo potential, (ii) all BHs (both the single and the binary) 
escaped the halo, and (iii) the single BH either escaped or 
remained wandering far out in the halo at the end of the 
run. The single escaped in ~15-20% of the runs in all cases. 
If we also count runs where it remained wandering through 
the halo for of order a Hubble time, this fraction increases 
to ~40%. Double escapes (of both the binary and the single 
BH) were very rare. We get more wandering BHs in set Dl, 
since a larger fraction of the released energy is apportioned 
to the escaper when it is relatively lighter, the dynamical 
friction time is longer, and the run time is shorter. Runs 
SA produced less wandering BHs since the binary pair more 
often coalesced before the intruder had a chance to harvest 
much of its energy. Wandering was also more common in 
set HI, due to the scaling discussed in the previous 

paragraph. The escape fraction of course depends on the 
depth of the galactic potential well. Given the uncertainty 
and scatter in the m b h — M^aio relation and specificity of the 
prescription adopted, we must expect these numbers to vary 
somewhat in studies with different halo or stellar density 
models. 

The entries under "final state" tell whether, at z = 0, 
the galactic center hosts (i) a stalled BH binary, (ii) a single 
BH, or (iii) no BHs (neither the single nor the binary has yet 
returned to the center by dynamical friction). About 50-70% 
of the runs ended with a binary at the galactic center whose 



gravitational radiation time exceeded the time until z = 0. 
This includes cases where (a) the single was ejected to large 
distance and the binary settled to the center before it hard- 
ened enough to coalesce by gravitational radiation, (b) when 
the inner binary coalesced during a close encounter the outer 
binary coalescence time exceeded the remaining run time, or 
(c) the single and binary remnant both returned to the cen- 
ter after an ejection and formed a bound pair with a long 
gravitational radiation time. In most of the remaining cases 
(30-50%) the run ended with a single BH at the galactic 
center, or a binary bound to coalesce before to- This oc- 
curred when (a) the single was ejected to large distance and 
the binary (or remnant) settled to the center after having 
hardened to the point of coalescence through some combina- 
tion of repeated interactions with the third BH and stellar 
dissipation, or (b) a new binary with a short gravitational 
radiation time formed following return from an ejection or 
coalescence during a close encounter. In only a small frac- 
tion (1-2%) of cases the run ended with the center empty of 
BHs. Note also that this happened most often in runs where 
the last merger occurred recently, so the total time spent 
with the center empty of BHs is smaller still. 

The next two entries give the mean and standard de- 
viation of the core "mass deficit" scoured out by the triple 
system, in units of the total BH mass rribh- For a galaxy 
modelled as an 77- model with stellar mass parameter M B , 
bulge scale radius a a and a break to inner slope 7 at r b , wc 
define the mass deficit Mdef by 



M, 



def 



4-k 



pn{r)r 2 dr - 



PT 

Jo 



(- 



-r - (- 



p(r)r dr 



■ 7 



rt + DM., (37) 



r + a a r b + a s 

where pb s is the stellar density at r b and D.M. denotes 
the corresponding dark matter terms. The mass deficits 
are highly scattered within each set of runs, with typical 
Mdef /m b h ~ 1.4 ± 0.5. More extreme mass ratios (runs 
Dl) tended to produce smaller cores, while a narrower mass 
range (runs D5) gave somewhat larger ones. The fraction of 
runs ending with very high mass deficits varied strongly with 
A m ; for instance 17% of cases ended with Mdef /rribh > 2 
for A m = 0.5, vs. only (11%, 4.4%) for A m =(0.3, 0.1). 
The large cores in set SA arose mostly from enhanced core 
scouring during the creation of the initial hard binary, and 
so are more a consequence of the ICs than of the triple in- 
teractions themselves. This sensitivity of Mdef to the binary 
stalling radius is an interesting point in its own right. The 
larger cores in runs HI probably arise from the higher mean 
number of ejections and smaller fraction of runs ending in 
immediate coalescence as the BH mass is decreased. 21% of 
the runs in this set ended with Mdef /rn b h > 2 and 8.7% 
ended with Mdef /rribh > 2.5. The subject of core scouring 
will be discussed in further detail in §3.5. 

Both the core scouring effect and the coalescence rate 
induced by the encounters are significantly reduced for the 
extreme mass ratios in set Dl. One must keep in mind, how- 
ever, that halo mergers with these mass ratios are much more 
frequent than those above A m = 0.3 (see Fig. 3), so the 
cumulative effect of these events may be as high or higher 
than that of encounters with near-equal masses. To quantify 
this statement our simulations would need to be embedded 
in a merger tree that follows the formation of triple systems. 
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Finally, the last five lines in Table 1 give the statistics 
of the condition which formally terminated the run: (i) the 
single escaped the halo and the binary (or remnant) set- 
tled to rest at the galaxy center; (ii) the escaper and binary 
remnant formed a stalled binary or coalesced (in a few per- 
cent of these cases a bound binary never actually formed; 
the pair coalesced suddenly upon a very close periapsis pas- 
sage from an unbound orbit); (iii) the maximum physical 
time tmax = to — tfo was reached (in these cases one or 
more BHs were left wandering through the halo at the end 
of the run); (iv) the maximum time for a close encounter 
(t c hn = 3 x 10 7 yrs for our canonical runs) was exceeded; 
or (v) the timestep went to zero or a limit on the number 
of timesteps was reached at some stage of the integration, 
which always occurred in <1% of cases. Runs terminating 
on condition (iv) or (v) were left out when computing the 
upper entries in the table. 

In this slew of runs we have varied only a few of the rele- 
vant parameters; one might also try, for instance, varying or 
adding scatter to the halo mass prescription, further steep- 
ening the stellar bulge profiles or adding a disk component, 
and exploring vastly different BH mass scales, in particular 
the much lower (~ 10 4-5 Mq) masses that may be relevant 
at high redshift. One of the advantages of our method is 
the relative ease of varying the model and ICs. This paper 
should be viewed as a work in progress, in which we have 
developed a method that can be applied to 3-BH systems 
in whatever context they may arise. Given the qualitative 
similarity of the outcomes in the runs we've performed so 
far, we will focus on the canonical (CN) runs in the more 
detailed presentation of our results. 

3.2 Efficient binary coalescence 

The inner binary begins at auard, where the gravitational 
radiation time is r gw ~ 10 13 ~ 15 yrs, in our canonical runs. It 
must shrink by a factor of ~10 before gravitational radiation 
can cause coalescence in a Hubble time, or by a factor of 
~40 for r gw to become comparable to the dynamical friction 
time. The intruder helps to bridge this gap in several ways: 
(a) direct hardening of the binary through repeated 3-body 
interactions, (b) enhanced stellar hardening by scattering 
of stars into the loss cone and motion of the binary about 
the nucleus, and (c) enhanced gravitational radiation losses 
due to thermalization of the eccentricity during the chaotic 
encounters and eccentricity growth via the Kozai resonance. 

We find that the combination of these mechanisms leads 
to coalescence of the inner binary within the time to — tfo 
between the merger that formed the triple system and z = 
in a large majority of the runs. It is intructive to distin- 
guish the systems that coalesce by "collision" during a close 
3-body encounter from those that gradually harden enough 
to coalesce within the time to — tfo, through the cumulative 
effect of repeated encounters and loss-cone refilling while the 
binary wanders about the nucleus. In runs CN, 23% of the 
systems undergo collision during the first encounter, and an- 
other 19% coalesce during later close encounters, for a net 
42% collision rate. Thus about half of the total coalescence 
efficiency arises from collisions during close encounters, and 
the other half comes from gradual hardening over the course 
of the simulation. Kozai oscillations account for most of the 
collisions during the first encounter, while in later encoun- 
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Figure 8. Distribution of coalescence times for the inner binary 
(upper panel) and new binary (lower panel), formed by the third 
BH and coalesced binary remnant. The lower plot includes only 
runs where a new binary formed, and not e.g. cases where the 
single escaped or was left wandering far from the galactic center 
at the end of the run. 

ters random eccentricity variations are more likely to result 
in coalescence, since the binaries are harder. Our numbers 
are reasonable based on anal ytic estimates of the collision 
rate in chaotic encounters (e.g. IValtonen fc Karttunenll2006l . 
chapter 11). 

Fig. 8 shows the distribution of binary coalescence 
times. The upper panel is for the inner binary, while the 
lower panel is for the new system formed by the third BH 
and binary remnant. In cases where coalescence occurred 
during the run, we plot the coalescence time recorded by 
the code. In other cases we plot t run + t gr ^ en d, where t run 
is the total r un time and t ar ,end is the time obtained by in- 
tegrating the iPetersl (1 19641 ) equations from the state at the 
end of the run to coalescence. The lower plot includes only 
those runs where the third BH ended up bound to the bi- 
nary remnant, excluding, for instance, cases where the single 
escaped the galaxy. In ~15% of the runs the new binary also 
coalesced within the time to — tfo. 

Under circumstances where the gap-crossing mecha- 
nisms discussed in §1.1 fail, the efficient coalescence in mas- 
sive triple systems provides a "last resort" solution to the 
final parsec problem. 

3.3 The 3-body interactions 

Though the close encounters take up only a small fraction 
of the physical time in our runs, it is the energy exchanges 
during these encounters that determine the large-scale BH 
dynamics. We now take a closer look at the 3-body dynamics 
in a few representative cases. 

In ~20% of the runs the binary swiftly coalesces dur- 
ing the first encounter, usually with the help of the Kozai 
resonance. Two examples of this are shown in Fig. 9. The 
time evolution of the inner and outer binary separations is 
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Figure 9. Two examples of rapid coalescence by Kozai oscilla- 
tions. Upper panels: Time evolution of the inner (red, lower) and 
outer (blue, higher) binary separations. Lower panels: Total grav- 
itational radiation power, averaged over Bulirsch-Stoer timesteps. 
mi and m2 are the masses of the binary members and mz is the 
mass of the intruder. 



plotted for two different runs in the upper panels. For a cir- 
cular orbit the separation would be roughly constant over an 
orbital period, or just a horizontal line in the figure. On the 
left the inner binary undergoes many Kozai oscillation cycles 
before coalescing. Observe that at the second-to-last eccen- 
tricity maximum, though it does not coalesce, the binary 
radiates away a large amount of energy and passes through 
the next eccentricity minimum with a significantly reduced 
semi-major axis. In the example on the right, the binary co- 
alesces after just one full Kozai cycle. The lower panels show 
the time evolution of the total gravitational radiation power, 
averaged over the Bulirsch-Stoer timesteps. Since the system 
starts on the verge of chaotic interactions where the outer 
to inner binary semi-major axis ratio is small (so that the 
quadrupolar approximation breaks down), we get "messy" 
Kozai oscillations which can give way to catastrophic eccen- 
tricity growth at an unpredictable time. 

Fig. 10 shows two examples of more complex runs. The 
left panels summarize the entire run, including all of the 
close encounters and ejections in between. Each call to the 
Chain code or the unperturbed binary integrator is sepa- 
rated by dashed vertical lines. The total time in each stage 
is normalized to unity in order to see the full history of the 
run at once, and not just distant ejections. The numbers on 
the plot are the actual times (in yrs) spent in each stage. 

Each run begins with a short period of secular evolution 
(illustrating the remarkable stability of hierarchical triples 
even slightly within the Mardling-Aarseth boundary). Dy- 
namical friction brings the intruder in a bit further to get 
chaotic interactions underway. This can be seen more clearly 
in the upper right panel, where we zoom in on the first close 
encounter at left. In this panel we also color-code the lines 
according to which two BHs instantaneously form the clos- 
est pair, to show the numerous exchanges that occur dur- 
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Figure 10. Two examples of longer runs. Left panels: Entire run, 
with time spent in each code stage normalized to unity. Actual 
times in years are indicated by the numbers on the plot. The red 
(inner binary) and blue (outer binary) portions show the close 
encounters, while the black portions show the calls to the RK4 
integrator. Upper right: Zoom in on the first close encounter in 
the run at left. Lines are color-coded according to which pair is 
closest, to highlight the exchanges. Lower right: Zoom in on the 
third close encounter of the lower run, showing also the total grav- 
itational radiation power averaged over Bulirsch-Stoer timesteps. 
See text for further explanation of this figure. 



ing close encounters. Large-amplitude Kozai oscillations are 
present in the first encounter of the lower run, but no oscilla- 
tions are seen in the upper run, where the initial inclination 
is below the critical angle. 

After the first encounter, in both runs the outer com- 
ponents suffer a few "near" (~0.1-1 kpc) ejections before 
they get shot out to kpc scales and come back by dynamical 
friction. In the upper run, the single goes out to ~10 kpc, 
then comes back and forms a bound pair with the former 
binary, which has coalesced in the meantime. The new bi- 
nary is highly eccentric (e ~ 0.9998) and quickly coalesces by 
gravitational radiation. In the lower run the single returns 
after the first kpc-scale ejection, strongly interacts with the 
binary one more time, and then escapes the galaxy. The bi- 
nary has a semi-major axis of 0.15 pc at the beginning of 
the final encounter, and its binding energy increases by 17% 
in the interaction, imparting a velocity of ~1400 km/s to 
the escaper. 

In the lower right two panels we focus on the third en- 
counter of this run, which was selected because a significant 
amount of energy was lost to gravitational radiation over its 
duration. We can see that the radiation loss occurred during 
two very close approaches, by two different BH pairs. If close 
3-body encounters between BHs are sufficiently common in 
the universe, such gravitational radiation spikes could be 
detectable with LISA. 

Fig. 11 shows the distribution of post-encounter veloci- 
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Figure 11. Distribution of ejection velocities. Upper panel: Single 
BH. Lower-panel: Binary COM. The total BH mass is typically 
~ 6 x 10 s M . 



ties, for the single and recoiling binary. Included in the plot 
are all close encounters in which (a) the binary and single 
are unbound at the end of the encounter; (b) the binding 
energy of the binary increases by at least 5% (to avoid nu- 
merous "glancing" encounters where 5F just barely exceeds 
the close encounter threshold), and (c) the encounter ends 
by SF falling below threshold (and not e.g. by coalescence of 
the binary or timing out). The dashed vertical line indicates 
the typical galactic (stellar bulge + halo) escape velocity, 
Vesc ~ 1400 km/s. We see that the single will sometimes es- 
cape the galaxy (or be ejected far out into the halo where 
the dynamical friction return time exceeds a Hubble time), 
but the binary will rarely go far. 

The upper panel of Fig. 12 shows the distribution of 
fractional changes in the binding energy of the binary dur- 
ing close encounters, 1 + A = 1 + (BEj - BE )/BE . The 
first encounter of each run is excluded from this plot, since 
it begins in a special hierarchical triple configuration and 
includes some dissipation by the dynamical friction used to 
bring in the intruder. The red line shows the best fit to the 
form /(l + A) = KA- 1/2 (l + A)~ 9/2 , with the normaliza- 
tion K depending on the mass rati o s and intruder velocity, 
predi cted by theory l|Heggiel Il975l ; IValtonen fc Karttunenl 
bOQtt ) . The lower panel shows the fraction of encounters with 
the relative energy radiated as gravitational waves in the en- 
counter greater than 1 + A gw — 1 + E gw / BEq. This shows 
that gravitational radiation plays a significant role in the dy- 
namics in only a few percent of the encounters ending in the 
escape of one component. Another ~ 20% of the encounters 
end in coalescence; gravitational radiation of course plays a 
significant role in all of these. 

Another point of interest is the statistics of the closest 
approach distances between two-body pairs during the en- 
counters. Besides their intrinsic significance, the distances 
of closest approach are related to the extent of tidal strip- 
ping of the BHs during the encounters. One can imagine 




Figure 12. Energy exchanges during close encounters. Upper 
panel: Distribution of the fractional change in the binding energy 
of the binary, 1 + A = 1 + (BEt — BEq) / BEq, in close encoun- 
ters. Blue points with 1.5cr Poisson error bars are the simulation 
data; the red line shows a comparison with the shape predicted 
by theory (see text). Lower panel: Fraction of encounters with the 
relative energy radiated as gravitational waves in the encounter 



greater than 1 + A s 



1 + E gw I BEq. Both panels include only 



encounters that end with the single BH unbound from the bi- 
nary, and exclude the first encounter of each run, which starts in 
a special hierarchical triple configuration and also includes some 
dissipation by the dynamical friction used to bring in the intruder. 



that if some stars, or even the inner portion of an accre- 
tion disk, remained bound to the individual BHs at the end 
of an encounter, then some ejected SMBHs might become 
observable. 

Since Bulirsch-Stoer timesteps are not at all infinites- 
imal (see §2.4), we cannot simply take the minimum over 
the discrete timesteps to be the closest approach distance. 
When the relative perturbation SF from the third body is 
small, one can obtain the periapsis distance analytically in 
the Keplerian two-body approximation. When 5F is larger, 
the minimum over the timesteps should give a better esti- 
mate since the timesteps tend to be smaller, but this state- 
ment is difficult to quantify. To construct the distance of 
closest approach in our simulations, we first identify any step 
where d\r\/dt = f ■ v switches sign from negative to positive 
and \r\ < 30000(r s i + r S 2) for any pair as a "passing step" 
containing a close approach. Here r = fi — r*2, v — v\ — V2, 
and r s i,2 are the Schwarzchild radii of the two pair members. 
We then iteratively bisect the timestep, evaluating f • v at 
each bisection to find the place where it switches sign until 
the distance between the two bodies converges to within a 
part in 10 6 . 

Fig. 13 shows the distribution of the tidal radius run 
at closest approach, for the closest and second closest pair. 
Since we are interested in observing ejected SMBHs, we only 
include encounters where the single escaped with a velocity 
above 940 km/s, the typical velocity needed to reach the 
stellar scale radius of ~3 kpc in our galactic model, rud is 
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Figure 13. Tidal radius r t id of the lighter pair member at closest 
approach. In this plot we include only encounters where one BH 
escaped at a speed above 940 km/s. Upper panels: r t id in units 
of the Schwarzchild radius of the lighter pair member, r s m i n . 
For reference, the red circles show the percent of the BH's mass 
contained within r t id in an ct-disk if the BH is accreting at the 
Eddington rate rh^ddy averaged over the encounters in each bin. 
Lower panels: r ti< i in pc. Red circles show the mass in stars (in 
Mq ) contained within r ti£ j in the Hernquist profile used to model 
the stellar component of the galaxy. Left panels: Closest pair; 
Right panels: Second closest pair. 



defined by the equation 

Grri2 



Sa t 



(d - r t id) 2 



Grrii 



(38) 



where mi is the reference mass being stripped (the smaller 
pair member), mi is the other point mass, and d is the dis- 
tance between mi and mi. We solve this polynomial equa- 
tion for r t id exactly rather than Taylor expanding about 
rud/d = to get the familiar expression r = (mi/2m,2) 1//3 
for the tidal radius, since rud/d is not generally small at 
closest approach for the near-equal mass problem at hand. 
The upper panels show r t id in units of the Schwarzchild ra- 
dius of the smaller BH. The red circles indicate the per- 
cent of this BH's mass contained within rud in an a- 
disk ilShakura fe Sunvaevlll973l ; lFrank. King, fe Rainell2002l; 
iNaravarJ 20031 ) accreting at the Eddington rate rhEdd, as- 
suming a — 0.1 and a radiative efficiency of e — 0.1. The 
lower panels give rud in pc, and here the red circles show the 
mass in stars within r t id in the Hernquist model representing 
the stellar bulge. 

The a-disk model assumes that the disk is not self- 
gravitating and breaks down as m — > rriEdd, so the red circles 
in the upper panels are merely to give the reader an idea of 
the bound mass scales associated with the approaches. The 
tidal approximation is a pessimistic estimate of the extent of 
the strippi ng since swift, one-time c lose passages would be 
impulsive (Bi nnev fe Trem ainc 1987j, chapter 7). We record 
only the single closest approach, so we cannot distinguish 
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Figure 14. Total time spent by the single BH (upper panel) and 
binary/remnant (lower panel) at various distances from the galac- 
tic center, averaged over all 1005 runs. The lowest bin includes all 
distances below 17 pc and the highest bin includes all distances 
above 10 5 pc. 



between such impulsive events and approaches that are part 
of periodic patterns in the trajectories. 

In a significant fraction of cases r t id J>10 4 r s encloses a 
substantial fraction of the BH's mass in accreting gas, so 
near-Eddi ngton ac c retion could continue for a duration of 
order the Salpetei (1964) time after the slingshot ejection 
l|Hoffman fe LoebT 20061 ). The enclosed stellar mass shown 
in the lower panels is never nearly comparable to the BH 
mass, but in most cases the escaper would drag some stars. 
In principle one can imagine one of these stars entering a 
giant phase and overflowing its Roche lobe, producing de- 
tectable accretion onto t he SMBH long after i ts ejection from 
the g alactic center (e.g. lHopman et al.| [2004; Ku ranov et al.l 
120071) . 



3.4 Distant Evolution and Binary Re-formation 

Slingshot ejections in triple encounters produce a popula- 
tion of "wanderin g" SMBHs in the halos of galaxies and in- 
terga lactic space (|Volonteri et al.ll2003al ; IVolonteri fe Pernal 
120051 ). Fig. 14 shows the total time spent by the single BH 
(upper panel) and binary /remnant (lower panel) at various 
distances from the galactic center, averaged over all 1005 
runs. The time spent in each distance bin is summed over 
the duration of each run, and if the system reaches a steady 
state at time t en d < to — i/o then the state of the system af- 
ter the run (until z = 0) is included. If a component escapes, 
then the time to — t en d is added to the highest bin; if it settles 
to the center then this time is added to the lowest bin. The 
single is found wandering at large distances nearly half the 
time, while the binary spends the vast majority of its time 
at the galactic center. Over all of our runs, the total fraction 
of the time spent with no SMBHs within 50 pc of the center 
since the formation of the halo hosting the triple system is 
only ~ 1%. Hence we expect the ejections in triple-BH en- 
counters at low redshift to produce very few nuclei empty 
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Figure 15. Distribution of parameters of "final state" binaries. 
The plot includes binaries that settled to the center after a single 
escape, and new binaries that formed from the single BH and in- 
ner binary remnant. Upper panel: Semi-major axis. Lower ; 
Eccentricity. 



of SMBHs. A cD galaxy cluster, having hosted several dry 
mergers, might contain up to a few naked SMBHs wandering 
through the cluster halo as a result of single ejections. 

The escaper remains wandering through the halo in only 
~40% of the runs. In the other cases dynamical friction 
brings it back to the center, where it becomes bound to 
the binary remnant and forms a new, hard binary. Fig. 15 
shows the semi-major axis (upper panel) and eccentricity 
(lower panel) distributions of the "final state" binaries in 
our simulations. This plot includes binaries formed when a 
pair coalesces during a close encounter and is replaced by a 
single BH with its COM coordinates; cases where the orig- 
inal binary never coalesces, but rather settles to the center 
and stalls after the single escapes; and binaries that form 
from the third BH and coalesced remnant after ejections. 

Whereas in the absence of triple encounters we would 
expect most SMBHs to sit around ahard, the encounters in- 
troduce a second population of stalled binaries at smaller 
separations. The eccentricities of the final binaries span the 
whole range from to 1. Note the peak at very high eccen- 
tricity, arising mostly from cases where the escaper rejoins 
the binary remnant from a radial orbit following a distant 
ejection, as in the run shown in the upper panel of Fig. 10. 
Many of the binaries in this peak are expected to coalesce 
quickly by gravitational radiation. This result has impor- 
tance for LISA if 3-body ejections are common enough, since 
the gravitational radiation signature of a highly eccentric bi- 
nary is quite different from that of a circular binary. An ec- 
centric binary radiates at all integer harmonics of the orbital 
frequency, so its spectral energy distribution peaks at higher 
frequencies, possibl y enabling the detection of higher-mass 
SMB H binaries (e.g. iPierro et ai1l200ll ; lEnoki fc Nagashimal 
l200tj ). 

However we caution the reader that the high- 
eccentricity coalescence rate appears to be sensitive to 
the dynamical friction and core updating prescriptions. 



SMBH binary eccentricity evolution is a delicate question 
to which si mulators have o btained widely discrepant an- 
swers (e.g. lAarseth I 120031: i Milosavlievic fc MerrittJ l200ll : 



iMerritt fc Mirosavlieviall2005h . While the conclusion that 
high-eccentricity population forms through distant ejections 
is robust, the question of whether these systems tend to 
coalesce or remain as an observable population of stalled 
high-eccentricity binaries may depend on the details. 



3.5 Core creation 

A number of studies have addressed the mark left on a 
stellar core by the hardenin g of one or more BH b i naries 
in independent succession ([Milosavlicvic & Me rrittl l200ll: 
Ravindranath et all 120021 ; IVolonteri et all l2003bl : iMerrittJ 
20061 ). To estimate the damage, consider a succession 
of mergers (M ,Mq) -> Mi, (Mi,M{) -> M 2 , 
(Mjv-i, M^r_!) — > Mjv, between galaxies containing BHs 
of mass (mo, mo), ... , (mjv-i, m' N _ 1 ), and having insuffi- 
cient gas for significant stellar cusp regeneration. Suppose 
that following each merger the BHs spiral in to ahard by dy- 
namical friction on the stars, then cross the gap from a na rd 
to a gw by some non-stellar mechanism, e.g. interaction with 
a modest amount of gas that ends up in the nucleus through 
tidal torques associated with the merger. The total energy 
deposited in the stellar core is roughly 
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where a nar d,i is the hardening radius of the BH binary 
formed in each merger, and the radius of influence o» n /,i 
is the radius containing abou t twice the m ass of the larger 
binary member in stars (e.g. IMerrittJ 120061 ) . Note that the 
right-hand side of (|39[) is dominated by the first term in the 
parentheses, so the precise definition of a in fj is not impor- 
tant. If the inner density profile of a galaxy flattens from 
dlnp/dlnr = 3 — rj to a shallower slope 7 within some core 
radius rt>, then we can define a core "energy deficit" by 
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(40) 



the difference between the binding energy of the galaxy 
with the density break and that of the same galaxy, but 
with the density profile outside the core extrapolated in- 
ward to the center. In this equation p — p s tars + phaio and 
& — $jiars + $haio denote the sums of the contributions to 
the density and gravitational potential from the stellar and 
dark matter halo components. The cross terms pstars^haio 
and phaio^ stars contribute about 10-20% of the total binding 
energy, while the halo-halo terms ar e negligible. We denote 
the outer slope by 3 — rj to match the lTremaine et"al1 (|l994h 
parameterization used in our galactic models. We can esti- 
mate the size of the core created by the cumulative scour- 
ing action of the BH binaries formed in the succession of 
mergers by equating Udef of equation (|40[) to Edep given by 
equation (|39[1 . 

The Udef = Edep prescription was introduced in order 
to estimate the extent of cusp destruction before the binary 
hardens. If stalling were prevented by sufficient scattering 
of stars into the loss cone, an analogous energy argument 
would grossly overestimate the size of the core scoured out 
as the binary decayed from a^ard to the separation where 



Triple black hole systems 19 



gravitational radiation could take over. This is because a 
hard binary loses energy by ejecting stars at high veloci- 
ties, often exceeding the escape velocity of the entire galaxy. 
Most of the energy released by the binary goes into excess 
kinetic energy of these hyper-velocity stars rather than heat- 
ing of the local medium. Equations (|39|l and (|40p capture the 
essence of the core scouring in the limit of weak encounters 
(dynamical friction), but for hard binaries we must view the 
cusp destruction as mass ejection rath er than energy injec- 
tion, o nce again following the work of iMikkola fc Valtonenl 
(1992) and lQuinlanl |l996l ). A hard SMBH binary is defined 
by the fact that it hardens at a constant rate, dE/dt — const. 
In the limit of very high orbital velocity (w 3> a) this im- 
plies that a constant mass in stars, comparable to the total 
BH mass, is ejected from the galactic center per e-folding of 
the binary semi-major axis, 



dM e 



dM e 



J ^ 0.5, 



(41) 



m bin dln(l/a) m bin d\n(E 3 
where Est is the energy transferred from the BH system to 
the stars. We can estimate the core damage due to mass 
ejection by equating the total mass ejected by the binary to 
Mdef as defined in equation (I37|l . 

Now suppose that instead of coalescing without further 
damaging the stellar core, the binary formed in the first 
merger in our sequence stalls at ahard until a third BH sinks 
in following the second merger. On the one hand, some en- 
ergy that would have been injected into the stars as the 
outer binary hardened may now instead be carried off as 
gravitational radiation or kinetic energy of a fast escaping 
BH, causing less damage to the stellar core than the decay 
of two separate binaries. On the other hand, the intruder 
may continue scattering stars into the loss cone well after 
the inner binary reaches ahard, and ejected BHs heat the 
core by dynamical friction as their orbits pass repeat edly 
through the dense nucleus dBovlan-Kolchin et al.ll2004l ). 

To quantify these considerations our code evolves the 
core radius r b and slope 7 along with the BH orbits to 
roughly account for the core heating and mass ejection 
caused by the triple systems. At the beginning of each run 
we initialize the core by injecting an energy 

_Gm\m2 G(mi +m2)m3 
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into the parent 77-model according to equation (|40[) . In runs 
where the inner binary starts at a > ahard we replace ahard,i 
with ai n u t i in equation (|42p . In runs where it starts at a < 
a-hard we also eject a mass 0.5mM„ ln(ahard,i/a init ,i) accord- 
ing to equation (|37[) before the start of the run. 

There is an ambiguity in the way we update the pro- 
file since energy may be injected (or mass may be ejected) 
either by increasing r b to make the core larger, or by de- 
creasing 7 to make it shallower. We resolved this ambiguity 
b y performing a rough fit to the 7 vs. y s Mdef /m b h data 
in iMerrita ([2006), to obtain the relation 

7 « -0.0281y 3 + 0.2451y 2 - 0.7094y + 1.000 (43) 

for r\ = 2, which gives sensible slopes for all y <^5. This 
relation at least has the desired property that 7 — > 3 — 77 as 
Mdef — > 0, but the slope becomes quite shallow toward large 



Mdef ■ The mass deficits are not sensitive to our prescription 
for 7. 

During the unperturbed binary integration, we incre- 
ment the energy injected at each timestep t — > t + At by 

2 pt+At 

AE in j = (Fdf,i • Vr)0{r cor e - r t )dt, (44) 

i=l Jt 

the work done on the single BH and binary COM by dynam- 
ical friction while the respective masses are located within 
a distance r core = 1.5rj of the galactic center. When the bi- 
nary is located within r core and has not yet coalesced or set- 
tled to the center and stalled, we also increment the ejected 
mass by 



AM ej = 0.5m Wn ln ^ ~ + J^ Est 



(45) 



at each timestep, where Eq is the binding energy of the 
binary at the beginning of the timestep and AE st is the 
change in binding energy due to stellar hardening during 
that step. We cannot simply use the semi-major axis incre- 
ment in equation (|45[1 since the binary may also have hard- 
ened by gravitational radiation during the timestep. At the 
beginning of each run, we also include the energy injected 
by dynamical friction as the intruder spirals in before the 
onset of chaotic interactions. 

Each time the total energy injected reaches 1% of Ehard 
or the total mass injected reached 1% of the total BH mass 
we update the core accordingly. At the end of each run we 
record the final Mdef , r b , and 7. Though both mass ejection 
and energy injection enter into our core growth algorithm, 
at the end both translate into a single effective Mdef as 
given by equation (|37[1 , which we record for comparison with 
observed galaxies. 

Fig. 16 compares our calculated mass deficits to 14 
cored, luminous elliptical galaxies with BH masses ranging 
from ~ 10 s — 3 x 10 9 Mq, with measured mass deficits. 
The upper panel shows data from our simulations. The blue 
(left) histogram is the distribution at the beginning of the 
runs, reflecting the heating of the core by dynamical friction 
on the BHs as they sink into the initial configuration with 
the inner binary at ahard, equation (142[) . This distribution 
also approximately represents the core damage expected for 
a single merger in which stellar hardening ceases at ahard- 
Note that a significant core {Mdef /mbh ~ 0.5) is scoured out 
even before the binary hardens and begins ejecting stars. 
The middle (black dashed) histogram shows the core pre- 
dicted for a series of two dry mergers, in both of which stellar 
hardening stops at ahard- The red (right) histogram is the 
distribution of cores at the end of our runs, reflecting the 
core scouring effect of the 3-body interac tions. Th e mass 
defi cits in the lower pane l are obtained in iGraharr] (I2004T ) 
and lFerrarese et alj (|2006l ) by fitting th e outer nuclear den- 
sity profile to a Sersic law (|Sersid ll968). and then subtract- 
ing a power-law fit to the inner core from the inward extrap- 
olation of the Sersic profile. The red dashed histogram in the 
lower panel shows the observed Md ef /mbh with m b h com - 
puted fr om the m b h — o~ rel ation of iTremaine et al. I (|2002h - 
However lLauer et al. I (|2006h argue that luminosity may be 
a better predictor of BH mass than a for the most lu- 
minous elliptical galaxies (Mv <, — 22), since their recent 
merger histories consisted mostly of passive (dissipation- 
less) mergers in which both the BH mass and luminosity 
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Figure 16. Comparison of calculated and observed mass deficits. 
Upper panel: Calculated mass deficits in units of the total BH 
mass, mi + m2 +JTI3. Blue (left): Beginning of simulation, based 
on core heating during inspiral to the initial configuration. Black 
dashed (middle): Net effect of two dry mergers in which core 
scouring stops at a^ ar( j. Red (right): End of simulation, reflect- 
ing net energy injection and mass ejection caused by close triple 
encounters. L ower panel: Ob s erved mass deficits in units of the 
SMBH mass llGrahaml 12004 iFerrarese et al.ll2006h . In the red 
dashed histogra m we det e rmine the BH masses from the m — cr re- 
lation of lTremaine et ah! J2002T) . for all entries. In th e black solid 
histog ram we instead use the m — My relation of lLauer et al.l 
(2006) for those galaxies with My < —22, and the dynamically 
measured BH masses in the four cases where they are available. 



Table 2. Mass deficits in gala xies with dynamicall y mea- 
sured SMBH mas ses. Refe rences: iGebhardt et al I 120031 (G03), 
Bower et all Il998l (B98), iMacchetto et al.l Il997l (M97), and 
Harms et ai]|l994l (H941 . 



Galaxy 


m b h/M Q 


M v 


M de f/m bh 


Reference 


NGC 4291 


3.1 x 10 8 


-20.64 


1.8 


G03 


NGC 4374 


1.6 x 10 9 


-22.28 


1.1 


B98 


NGC 4486 


3.0 x 10 9 


-22.71 


2.9 


M97, H94 


NGC 4649 


2.0 x 10 9 


-22.51 


1.1 


G03 



are simply additive. The rribh — o relation is thought to 
arise from self-r egulation of accretion onto the SMBH i n 
gaseous mergers i|Silk fc Reel 19981 ; IWvithe fc Loebl|2003bh 
which does not apply in this context. lLauer et alj ~ 20061 ) 
also show that an extrapolation of the rribh — L relation to 
the highest luminosities is more consistent with the observed 
r C ore — rnth relation and provides a better match between 
the 2 = SMBH space density and the quasar population 
seen at z ~ 2 for reasonable quasar duty cycles. In the black 
histogram, we used the observed BH mas ses for the four 
cases with dynamical mass measurements l|Gebhardt et al 



2003) ; iBower et all 19981 ; IMacchetto et al.lll997l ; lHarms et al 
994 ). For the rest of the galaxies we used the lLauer et al 
2006) rribh — My relation to esti mate the BH ma s ses in 

galaxies with M v < -22, and the iTremaine et all (|2002l ) 



rribh — relation for those with Mv > —22. Showing both 
plots gives an idea of how much the mass deficits vary with 
rribh estimator. 

In set CN of the runs resulted in cores with 

Mdef /rribh > 2. In some rare runs where the binary was 
ejected to a large distance and then brought back by dy- 
namical friction, we obtained even higher mass deficits (up 
to Mdef /rribh = 3 — 4). The binary is more efficient at core 
scouring than the single BH since it is more massive. If 
each independent binary inspiral adds ~ 0.5rribh to the mass 
deficit, then the mean Mdef enhancement of ~ 0.5iribh due 
to the triple encounters is equivalent to one extra merger in 
the system's history. An Mdef one standard deviation above 
the mean is equivalent to two extra mergers. 



4 DISCUSSION AND CONCLUSIONS 

Triple-SMBH systems in galactic nuclei produce a range of 
phenomena and signatures rather different from those ex- 
pected if no more than two SMBHs occupy them at a time. 
We have developed an efficient numerical method for follow- 
ing the evolution of 3-body systems in the centers of galaxies, 
and used it to explore the outcomes of such encounters in 
massive elliptical galaxies at low redshift. 

We find a high efficiency of SMBH coalescence due to 
the encounters, providing a "last resort" solution to the fi- 
nal parsec problem. There is, however, a caveat in extending 
this result immediately to all BH masses. If we define a esc 
to be the binary semi-major axis where escape of one BH 
first becomes likely (Grribi„/ [3a eS c = v^ sc where f3 is a factor 
of order 10 for a Hernquist profile), then since Ve 3C oc a we 
have ciese oc nibin/a 2 oc rrCd* if mu n obeys the m b h ~ a re- 



lation rribv 



c/a g 



-1/4 



In other words, at 



A,g W UV. Hl-bi n 

smaller BH masses, the lightest BH is more likely to escape 
the galaxy before driving the binary to coalescence by grav- 
itational radiation. By focusing on massive galaxies we have 
chosen the systems where the binary is least likely to coa- 
lesce by other means (e.g. gas or massive perturbers), and 
most likely to coalesce in the next merger with the help of 
3-BH interactions. We may address the efficiency of triple- 
induced coalescence in much smaller-mass systems in future 
studies. 

We find that close triple encounters can produce a 
population of high-eccentricity binaries, whose gravitational 
radiation signal could potentially be observable by LISA. 
Such signals originate from Kozai oscillations in hierarchi- 
cal triples at high initial inclinations and highly eccentric 
binaries formed following distant ejections. As the eccentric- 
ity increases, the radiation spectrum peaks at progressively 
higher harmonics of the fundamental fr equency, approach- 
ing a nearly flat spectru m as e — > 1 (|Pierro et al.l l200ll ; 
lEnoki fc Nagashimall2fJ06l ). A circular 10 8 " 9 M Q BH binary 
remains below the band of frequencies (~ 10 -4 — 10 _1 Hz) 
detectable by LISA throughout its inspiral, but the occur- 
rence of high-eccentricity coalescences could extend LISA's 
sensitivity into this mass range, or lengthen the duration of 
its sensitivity to ~ 10 6-7 Mq events. A highly eccentric bi- 
nary produces a "spiky" waveform that looks q uite different 
from that of a circular system (see Fig. 7 in iPierro et al.l 
2001). Gravitational radiation "spikes" at very close ap- 
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proaches during chaotic 3-body interactions could also pro- 
duce radiation bursts detectable by LISA. 

If triple encounters are indeed limited to massive sys- 
tems at low redshift, then the importance of these considera- 
tions is limited by the expected event rate in this mass range, 
assuming efficient coale scence. This rate is highly uncertain, 
ranging from ~ l/ vr jSesana et al.l 120051 ) to ~ l/1000yrs 
l|Rhook fc Wvithell2005l ) depending on the merger and BH 
population model adopted. If 3-BH systems occur in other 
contexts, e.g. IMBHs in galactic nuclei or star clusters, then 
the phenomena we have discussed may be observationally 
relevant even if the high-mass SMBH event rate is low. A 
detailed look at the gravitational waveforms expected from 
3-body encounters and their expected detection rates is an 
interesting topic for a future study. 

The slingshot ejections in triple encounters produce a 
population of "wandering" SMBHs in and outside the halos 
of galaxies. In systems that have undergone several major 
dry mergers (e.g. cD galaxy systems), one might expect a 
few such ejected SMBHs to be floating in the vicinity. As of 
yet, no probable way of observing these wandering BHs has 
been proposed!]]. In principle one can imagine a star bound to 
the ejected SMBH entering a giant phase and overflowing its 
Roche lobe, producing some accretion onto the SMBH and 
an observable flare. Single ejections could also in principle 
affect BH-bulge correlations such as the m^-ff relation, but 
since it is the lightest BH that gets ejected this effect would 
fall well within the observed scatter in the correlations for 
just one or two ejection events. 

Triple interactions in galactic nuclei can have a large 
effect on the expected properties of stable SMBH binaries 
in the local universe. While many models of binary forma- 
tion predict mostly circular binaries around a^ard, 3-body 
encounters produce binaries at all eccentricities. They also 
create a population of stalled binaries at separations signif- 
icantly smaller than a^ard but still larger than a gm , as does 
any partial gap-crossing mechanism. 

Better measurements and statistics on the mass deficits 
in cored elliptical galaxies may provide clues on the history 
of the nuclear SMBH activity in these systems. Triple BH 
encounters produce a highly scattered distribution of core 
sizes, with mass deficits up to ~ 2x higher than expected for 
successive binary coalescences. The apparent peak at mass 
deficits of ~ 0.5 — 1 times the nuclear BH mass in observed 
cores may very tentatively hint that multiple-BH encounters 
are not the norm in these systems. This signature of binary 
or multiple-BH activity is appealling because (a) its duty 
cycle is the lifetime of the galaxy; (b) it is present whether 
binary pairs stall or coalesce; and (c) it can be observed 
even in the complete absence of radiative activity, such as 
disk accretion or jet production. However the interpreta- 
tion of galaxy cores is complicated by multiple mergers, the 
possibility of partial stellar cusp regeneration from traces 
of cold gas, and observational complications such as pro- 
jection effects in nonspherical galaxies and optimizing the 
fitting/extrapolation algorithm to best represent the mass 
deficit. There is a need for theoretical studies on the cores 
produced by SMBH mergers in triaxial galaxies, since triax- 



1 Gravitational lensing is difficult to search for without knowing 
in advance the location of the BHs. 



ility seems to be the most likely candidate for a gap-crossing 
mechanism in dry mergers between gas-poor, giant ellipti- 
cals. Inferring the nuclear histories of galaxies from their 
observed core properties will likely be a topic of much inter- 
est in the future. 
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